diff --git a/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibT0/MdtCalibT0/ADCMTHistos.h b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibT0/MdtCalibT0/ADCMTHistos.h
new file mode 100644
index 0000000000000000000000000000000000000000..b2ad4acb9fd0f20259fa20b0f557535abb19ed8d
--- /dev/null
+++ b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibT0/MdtCalibT0/ADCMTHistos.h
@@ -0,0 +1,120 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+// ADCMTHistos.h
+// Header file for class ADCMTHistos
+///////////////////////////////////////////////////////////////////
+// (c) ATLAS Detector software
+///////////////////////////////////////////////////////////////////
+//
+///////////////////////////////////////////////////////////////////
+
+#ifndef ADCMTHistos_H
+#define ADCMTHistos_H
+
+//c - c++
+#include <iostream>
+
+//root
+#include "TH1.h"
+/*#include "TF1.h"
+#include "TDirectory.h"*/
+//class TH1F;
+class TF1;
+class TDirectory;
+
+//this
+#include "T0MTSettings.h"
+
+
+namespace MuonCalib {
+
+/**@class ADCMTHistos
+ Histogram and fitter class for drift time and pulsehight spectra
+ The rising slope is fitted by a fermi-function: @f$f(t)=r_{u,0} + \frac{A}{1+exp((t_0 - t)/T_0}@f$.
+ The falling slope is fitted by @f$g(t)=r_{u,max} + \frac{e(t,a,b)}{1+exp((t_{max} - t)/T_{max}}@f$ where @f$e(t, a, b) = a e^{bt}@f$.
+ A pattern recognition determines the fit ranges and the parameters @f$r_u@f$, @f$A@f$, @f$a@f$ and @f$b@f$.
+ @author Felix.Rauscher@Physik.Uni-Muenchen.De
+ @date June 2006
+
+}
+
+*/
+class ADCMTHistos
+ {
+ public:
+//---------------------------constructor----------------------------------------
+ /** Default Constructor */
+ inline ADCMTHistos() :
+ m_adc(NULL),
+ m_id(-1),
+ m_adc_fit(NULL),
+ m_adc_ok(false)
+ {}
+ /** Initializing constructor
+ @param id tube id
+ @param settings t0-fit settings: settings will be asked about histogram binning
+
+ */
+ inline ADCMTHistos(int id, const T0MTSettings * settings, const char *hname=NULL):
+ m_adc(NULL),
+ m_id(-1),
+ m_adc_fit(NULL),
+ m_adc_ok(false)
+ {
+ Initialize(id, settings, hname);
+ }
+//---------------------------public member functions----------------------------
+ /** Initialize class
+ @param id tube id
+ @param settings t0-fit settings: settings will be asked about histogram binning
+
+ */
+ void Initialize(int id, const T0MTSettings * settings, const char * hname=NULL);
+
+ /** get adc spectrum */
+ inline TH1F * GetADCSpec() {return m_adc;}
+
+ /** fill adc value
+ @param a adc value
+
+ */
+ void FillA(double a);
+
+
+ /** return tube id
+
+ */
+ inline int Id() const {return m_id;}
+
+ /** fit adc*/
+ bool FitAdc();
+
+ /** returnd function fitted to adc-spectrum */
+ inline TF1 * GetAdcFunction() const
+ {
+ return m_adc->GetFunction("adc_fun");
+ }
+ inline bool AdcOk() const
+ {
+ return m_adc_ok;
+ }
+
+ private:
+//---------------------------private data members-------------------------------
+ //! pulse height spectrum
+ TH1F *m_adc;
+ //! tube id;
+ int m_id;
+ //! function which is fitted to the adc-spectrum
+ TF1 *m_adc_fit;
+ //! is true if the adc-fit is ok;
+ bool m_adc_ok;
+ //! Pointer to settings class
+ const T0MTSettings *m_settings;
+ };
+
+}
+
+#endif
diff --git a/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibT0/MdtCalibT0/HistogramId.h b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibT0/MdtCalibT0/HistogramId.h
new file mode 100644
index 0000000000000000000000000000000000000000..bddf667f14531441728b862603583bace1076464
--- /dev/null
+++ b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibT0/MdtCalibT0/HistogramId.h
@@ -0,0 +1,78 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#ifndef HistogramId_H
+#define HistogramId_H
+
+//c - c++
+#include <map>
+#include <string>
+#include <utility>
+
+namespace MuonCalib {
+
+class MuonFixedId;
+
+/** @class HistogramId
+Identifier class for drift time histograms. The class can be switched to the identification of single tubes, and different types of groupes of tubes.
+
+@author Felix.Rauscher@Physik.Uni-Muenchen.de
+*/
+
+class HistogramId
+ {
+ public:
+ /** constructor - does nothing*/
+ inline HistogramId():m_id(0,0) {}
+
+ /** initialize
+ @param id MuonCalibFixedId of the tube
+ @param sort_by (HistogramId :: TUBE, HistogramId :: LAYER ..)
+ */
+ void Initialize(const MuonFixedId &id, int sort_by);
+ /** comparision operator defined so that this class can
+ be used as a map index. The behaviour depends on the setting of the
+ sort_by parameter of the Initialize function.
+ */
+ inline bool operator > (const HistogramId & other) const
+ {
+ return m_id > other.m_id;
+ }
+ /** comparision operator defined so that this class can
+ be used as a map index. The behaviour depends on the setting of the
+ sort_by parameter of the Initialize function.
+ */
+ inline bool operator < (const HistogramId & other) const
+ {
+ return m_id < other.m_id;
+ }
+ /** valid values of the sort_by argument of the Initialize function */
+ static const int TUBE=4,
+ LAYER=2,
+ MULTILAYER=1,
+ CHAMBER=0,
+ MEZZ_CARD=3;
+ /** get numberic id */
+ inline int getIdInt() const
+ {
+ return m_id.second;
+ }
+ /** get ascii histogram name */
+ inline const std::string & HistogramName() const
+ {
+ return s_histogram_names[m_id];
+ }
+ private:
+ //! integer identity - value depends on the sort_by argument of the Initialize-Function
+ std::pair<int, int> m_id;
+ //! ascii histogram names
+ static std::map<std::pair<int, int>, std::string> s_histogram_names;
+
+ };
+
+
+}
+
+
+#endif
diff --git a/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibT0/MdtCalibT0/MTT0PatternRecognition.h b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibT0/MdtCalibT0/MTT0PatternRecognition.h
new file mode 100644
index 0000000000000000000000000000000000000000..dc41d0924574a4fdd78c18be2ce96ded0b6fe948
--- /dev/null
+++ b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibT0/MdtCalibT0/MTT0PatternRecognition.h
@@ -0,0 +1,122 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+// MTT0PattternRecognition.h
+// Header file for class VariableBinwidthHistogram
+///////////////////////////////////////////////////////////////////
+// (c) ATLAS Detector software
+///////////////////////////////////////////////////////////////////
+//
+///////////////////////////////////////////////////////////////////
+
+
+#ifndef MTT0PATTTERNRECOGNITION_H
+#define MTT0PATTTERNRECOGNITION_H
+
+//c - c++
+//#include "iostream"
+class TH1F;
+
+//root
+#include "TH1.h"
+
+//this
+#include "VariableBinwidthHistogram.h"
+#include "T0MTSettings.h"
+
+namespace MuonCalib {
+
+/** @class MTT0PattternRecognition
+ Performs the pattern recognition for the MT t0-fit routine. Fir ranges and start values are estimated
+ @author Felix.Rauscher@Physik.Uni-Muenchen.De
+ @date February 2006
+
+*/
+
+
+class MTT0PatternRecognition {
+
+ public:
+//-------------------------constructors-----------------------------------------
+ /** Default constructor*/
+ inline MTT0PatternRecognition(): m_settings(NULL){}
+
+//-------------------------public member functions------------------------------
+ /** Initialize class - returns true if pattern recognition was successfull
+ @param hist Histogram which is to be fitted
+
+ */
+ inline bool Initialize( TH1F *hist, const T0MTSettings *settings)
+ {
+ m_settings=settings->T0Settings();
+ m_draw_debug_graph=settings->DrawDebugGraphs();
+ m_error=!m_vbh.Initialize(hist, m_settings->VBHBinContent(), m_settings->MaxBinWidth());
+ double scale_min;
+ if(!m_error) scale_min=estimate_height(hist);
+ if(!m_error) estimate_background(hist, scale_min);
+ return !m_error;
+ }
+
+ /** get the background level */
+ inline double GetBackground() const {return m_background;}
+
+ /** get height */
+ inline double GetHeight() const {return m_height;}
+
+ /** get estimated t0*/
+ inline double GetEstimatedT0() const {return m_t0_est;}
+
+ /** get fit range */
+ inline double GetFitRangeMin() const {return m_fit_min;}
+
+ /** get fit range */
+ inline double GetFitRangeMax() const {return m_fit_max;}
+
+ /** return error flag */
+ inline bool GetError() const {return m_error;}
+
+ private:
+//------------------------private data members----------------------------------
+
+ //! settings
+ const T0MTSettingsT0 * m_settings;
+ bool m_draw_debug_graph;
+
+ //!background level
+ double m_background;
+
+ //!height
+ double m_height;
+
+ //!fit range
+ double m_fit_min, m_fit_max;
+
+ //!t0 estimate
+ double m_t0_est;
+
+ //!Variable binwidth histogram
+ VariableBinwidthHistogram m_vbh;
+
+ //! error flag
+ bool m_error;
+
+//-----------------------private member functions-------------------------------
+
+ /**estimates the background level
+ @param hist input histogram
+ @param scale_min lower end of the scale region as returned from estimate_height()
+
+ */
+ bool estimate_background( TH1F *hist, double scale_min);
+
+ /**estimates the height of the spectrum. It returns the lower end of the scale region
+ @param hist input histogram
+
+ */
+ double estimate_height( TH1F *hist);
+
+ };
+
+}
+#endif
diff --git a/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibT0/MdtCalibT0/MTTmaxPatternRecognition.h b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibT0/MdtCalibT0/MTTmaxPatternRecognition.h
new file mode 100644
index 0000000000000000000000000000000000000000..1a3a789eb5bf5672c5af37b4b691bc2cda3137c3
--- /dev/null
+++ b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibT0/MdtCalibT0/MTTmaxPatternRecognition.h
@@ -0,0 +1,120 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+// MTTmaxPattternRecognition.h
+// Header file for class VariableBinwidthHistogram
+///////////////////////////////////////////////////////////////////
+// (c) ATLAS Detector software
+///////////////////////////////////////////////////////////////////
+//
+///////////////////////////////////////////////////////////////////
+
+
+#ifndef MTTMAXPATTTERNRECOGNITION_H
+#define MTTMAXPATTTERNRECOGNITION_H
+
+//c - c++
+#include "iostream"
+
+//root
+//#include "TH1.h"
+class TH1F;
+
+//this
+#include "VariableBinwidthHistogram.h"
+#include "T0MTSettings.h"
+
+namespace MuonCalib {
+
+/** @class MTTmaxPattternRecognition
+ Performs the pattern recognition for the MT tmax-fit routine. Fir ranges and start values are estimated
+ @author Felix.Rauscher@Physik.Uni-Muenchen.De
+ @date February 2006
+
+*/
+
+
+class MTTmaxPatternRecognition {
+
+ public:
+//-------------------------constructors-----------------------------------------
+ /** Default constructor*/
+ inline MTTmaxPatternRecognition() :
+ m_settings(NULL),
+ m_error(false) {}
+
+//-------------------------public member functions------------------------------
+ /** Initialize class
+ @param hist Histogram which is to be fitted
+ @param t0 comes from the t0 fit
+
+ */
+ bool Initialize( TH1F *hist, double t0, const T0MTSettings * settings);
+
+ /** get the background level */
+ inline double GetBackground() const {return m_background;}
+
+ /** get parameter a in exp-function representing the end of the spectrum */
+ inline double GetA() const {return m_a;}
+
+ /** get parameter a in exp-function representing the end of the spectrum */
+ inline double GetB() const {return m_b;}
+
+ /** get estimated t0*/
+ inline double GetEstimatedTMax() const {return m_tmax_est;}
+
+ /** get fit range */
+ inline double GetFitRangeMin() const {return m_fit_min;}
+
+ /** get fit range */
+ inline double GetFitRangeMax() const {return m_fit_max;}
+
+ /** return error flag */
+ inline bool GetError() const {return m_error;}
+ private:
+//------------------------private data members----------------------------------
+
+ //! settings
+ const T0MTSettingsTMax *m_settings;
+ bool m_draw_debug_graph;
+
+ ///background level
+ double m_background;
+
+ ///parameters a, b in a*exp(b*x)
+ double m_a, m_b;
+
+ ///fit range
+ double m_fit_min, m_fit_max;
+
+ ///t0 estimate
+ double m_tmax_est;
+
+ ///Variable binwidth histogram
+ VariableBinwidthHistogram m_vbh;
+
+ ///t0
+ double m_t0;
+
+ /// error flag
+ bool m_error;
+
+//-----------------------private member functions-------------------------------
+
+ /**estimates the background level
+ @param hist input histogram
+
+ */
+ bool estimate_background( TH1F *hist);
+
+ /**estimates the height of the spectrum. The height is representet by an exponetial function a*exp(b*t). It returns the upper end of the height region.
+ @param hist input histogram
+
+ */
+ bool estimate_height( TH1F *hist);
+
+ };
+
+}
+#endif
diff --git a/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibT0/MdtCalibT0/ReflexHeaders.h b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibT0/MdtCalibT0/ReflexHeaders.h
new file mode 100644
index 0000000000000000000000000000000000000000..e61c9d4e9b3853867cc5d2057bce0865cdb90453
--- /dev/null
+++ b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibT0/MdtCalibT0/ReflexHeaders.h
@@ -0,0 +1,9 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#include "MdtCalibT0/T0MTHistos.h"
+#include "MdtCalibT0/T0MTSettings.h"
+#include "MdtCalibT0/T0MTSettingsT0.h"
+#include "MdtCalibT0/T0MTSettingsTMax.h"
+
diff --git a/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibT0/MdtCalibT0/T0CalibrationClassic.h b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibT0/MdtCalibT0/T0CalibrationClassic.h
new file mode 100644
index 0000000000000000000000000000000000000000..3cf194b56257286fe43f98884b53041c52a95f0f
--- /dev/null
+++ b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibT0/MdtCalibT0/T0CalibrationClassic.h
@@ -0,0 +1,176 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+// T0CalibrationClassic.h
+// Header file for class T0CalibrationClassic
+///////////////////////////////////////////////////////////////////
+// (c) ATLAS Detector software
+///////////////////////////////////////////////////////////////////
+//
+///////////////////////////////////////////////////////////////////
+
+#ifndef T0CALIBRATIONCLASSIC_H
+#define T0CALIBRATIONCLASSIC_H
+#include <iostream>
+#include <vector>
+
+#include "MdtCalibInterfaces/IMdtCalibration.h"
+#include "MdtCalibData/MdtTubeFullInfoContainer.h"
+#include "MdtCalibData/MdtTubeFitContainer.h"
+
+
+class TH1;
+class TH2;
+class TFile;
+class TDirectory;
+
+namespace MuonCalib {
+
+class MuonCalibSegment;
+class IMdtCalibrationOutput;
+
+/**@class T0ClassicSettings
+ Settings for the T0 calibration (histogram booking and fitting parameters)
+ @author domizia.orestano@cern.ch
+ @date June 2005
+
+*/
+class T0ClassicSettings{
+ public:
+ /** constructor */
+ T0ClassicSettings( double minAdc, double maxAdc, int adcBins,
+ double minTime, double maxTime, int timeBins,
+ bool fitTime , int minEntries, int initParam,
+ int numParams, double fitParams[], double chiMax, int printLevel) :
+ m_minAdc(minAdc),m_maxAdc(maxAdc),m_binAdc(adcBins),
+ m_minTime(minTime),m_maxTime(maxTime),m_binTime(timeBins),
+ m_fitTime(fitTime),m_entries(minEntries),
+ m_initParamFlag(initParam),m_numParams(numParams),
+ m_params(fitParams), m_chi2max(chiMax), m_printLevel(printLevel) { }
+ /** Access functions to values of private settings */
+ /** @return m_minAdc, minimum for Adc histogram */
+ double minAdc() const {return m_minAdc;};
+ /** @return m_maxAdc, maximum for Adc histogram */
+ double maxAdc() const {return m_maxAdc;};
+ /** @return m_binAdc, number of bins for Adc histogram */
+ int binAdc() const {return m_binAdc;};
+ /** @return m_minTime, minimum for Time histogram */
+ double minTime() const {return m_minTime;};
+ /** @return m_maxTime, maximum for Time histogram */
+ double maxTime() const {return m_maxTime;};
+ /** @return m_binTime, number of bins for Time histogram */
+ int binTime() const {return m_binTime;};
+ /** @return m_fitTime, flag activating the fitting of Time histogram */
+ bool fitTime() const {return m_fitTime;};
+ /** @return m_entries, the minimum number of entries to perform an histogram fit */
+ int entries() const {return m_entries;};
+ /** @return m_initParamFlag: false if intial parameters are fixed, true if
+ they have to be searched */
+ int initParamFlag() const {return m_initParamFlag;};
+ /** @return number of parameters to be fit */
+ int numParams() const {return m_numParams;};
+ /** @return the pointer to the vector of parameters */
+ double * params() const {return m_params;};
+ /** @return the maximum values of chi2 to accept the fit result */
+ double chi2max() const {return m_chi2max;};
+ /** @return the print level */
+ int printLevel() const {return m_printLevel;};
+ /** a method to dump all the settings */
+ void print() const {
+ std::cout <<"T0ClassicSettings: "<<std::endl
+ <<"booking of Adc histos: min max bins " <<minAdc()<<" "<<maxAdc()<<" "<<binAdc()<<std::endl
+ <<"booking of Time histos: min max bins " <<minTime()<<" "<<maxTime()<<" "<<binTime()<<std::endl
+ <<"fitting of Time histos: active minentries paramflag chi2cut nparams"
+ <<fitTime()<<" "<<entries()<<" "<<initParamFlag()<<" "<<chi2max()<<" "<<numParams()<<std::endl;
+ std::cout<<"Parameters: ";
+ for(int i=0;i<numParams();i++)std::cout<<m_params[i]<<" ";
+ std::cout<<std::endl;
+ };
+
+ private:
+ /** Private settings */
+ double m_minAdc; //!< Adc spectrum low edge
+ double m_maxAdc; //!< Adc spectrum high edge
+ int m_binAdc; //!< Adc spectrum number of bins
+ double m_minTime; //!< Time spectrum low edge
+ double m_maxTime; //!< Time spectrum high edge
+ int m_binTime; //!< Time spectrum number of bins
+ bool m_fitTime; //!< flag to select or deselect the fit to Time Spectra
+ int m_entries; //!< minimum number of entries per spectrum fitting
+ int m_initParamFlag; //!< fix or search initial parameters for spectrum fit
+ int m_numParams; //!< number of paramaters for spectrum fit
+ double * m_params; //!< initial parameters for spectrum fit
+ double m_chi2max; //!< normalized chi2 cut
+ int m_printLevel; //!< output level flag (1 for verbose)
+};
+
+/**@class T0ClassicHistos
+ Tube histograms used in T0 calibration.
+ @author domizia.orestano@cern.ch
+ @date June 2005
+*/
+class T0ClassicHistos {
+ public:
+ TH1* time; //!< time spectrum
+ TH1* adc; //!< adc spectrum
+ //TH2* adc_vs_time;
+ int id; //!< tube identifier
+};
+
+/**@class T0CalibrationClassic
+ Implementation of a T0 calibration using the classical approach.
+ @author domizia.orestano@cern.ch
+ @date June 2005
+*/
+class T0CalibrationClassic : public IMdtCalibration {
+
+ public:
+ /** constructor
+ * @param[in] name of the region/chamber to be calibrated
+ * @param[in] pointer to settings vector
+ */
+ T0CalibrationClassic( std::string name, const T0ClassicSettings* settings );
+ /** destructor */
+ ~T0CalibrationClassic();
+
+ bool handleSegment( MuonCalibSegment& seg ); //!< fill tube spectra
+
+ void setInput( const IMdtCalibrationOutput* input ); //!< unused
+
+ bool analyse(); //!< extract parameters from spectra
+
+ bool converged() const; //!< return m_converged (always false?)
+
+ /** @return the calibration results */
+ const IMdtCalibrationOutput* getResults() const;
+
+const IMdtCalibrationOutput* analyseSegments( const std::vector<MuonCalibSegment*> & segs );//!< new interface function
+
+ private:
+
+ T0ClassicHistos * bookHistos(unsigned int idtube); //!< booking of histograms
+
+ T0ClassicHistos * getHistos(unsigned int idtube); //!< retrieve pointer for tube idtube histograms
+
+ void doTimeFit(T0ClassicHistos*,MdtTubeFitContainer::SingleTubeFit&,MdtTubeFitContainer::SingleTubeCalib&); //!< fit time spectrum
+
+ void doAdcFit(T0ClassicHistos*, MdtTubeFitContainer::SingleTubeFit&,MdtTubeFitContainer::SingleTubeCalib&); //!< fit adc spectrum
+
+ void searchParams(TH1 * h, double * p, int np); //!< estimate initial parameters for time spectrum fit from the spectrum itself
+
+ const T0ClassicSettings* m_settings; //!< pointer to the settings
+ bool m_converged; //!< convergence status
+ std::string m_name; //!< calibration region name
+ int m_currentItnum; //!< current iteration (always 1?)
+ TFile * p_file; //!< pointer to the histogram file
+ TDirectory * m_regiondir; //!< pointer to the ROOT directory
+
+ std::vector<T0ClassicHistos*>p_histos; //!< vector of pointers tube histograms
+ MdtTubeFitContainer* m_result; //!< tube constants
+ bool m_delete_settings;
+};
+}
+
+#endif
+
diff --git a/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibT0/MdtCalibT0/T0CalibrationMT.h b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibT0/MdtCalibT0/T0CalibrationMT.h
new file mode 100644
index 0000000000000000000000000000000000000000..e9db0a4673a585947e95c594dc82fac298e59c7f
--- /dev/null
+++ b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibT0/MdtCalibT0/T0CalibrationMT.h
@@ -0,0 +1,118 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+// T0CalibrationMT.h
+// Header file for class T0CalibrationMT
+///////////////////////////////////////////////////////////////////
+// (c) ATLAS Detector software
+///////////////////////////////////////////////////////////////////
+//
+///////////////////////////////////////////////////////////////////
+
+#ifndef T0CALIBRATIONMT_H
+#define T0CALIBRATIONMT_H
+#include <iostream>
+#include <string>
+#include <vector>
+#include <map>
+#include <set>
+
+#include "MdtCalibInterfaces/IMdtCalibration.h"
+#include "MdtCalibData/MdtTubeFitContainer.h"
+#include "MdtCalibData/MdtTubeFitContainer.h"
+#include "MuonCalibStandAloneBase/NtupleStationId.h"
+//#include "T0MTSettings.h"
+//#include "T0MTHistos.h"
+//#include "ADCMTHistos.h"
+//#include "HistogramId.h"
+
+class TH1;
+class TFile;
+class TDirectory;
+
+namespace MuonCalib {
+
+class MuonCalibSegment;
+class IMdtCalibrationOutput;
+class MuonFixedId;
+class T0MTSettings;
+class T0MTHistos;
+class T0ADCHistos;
+class ADCMTHistos;
+class HistogramId;
+class MdtRelativeTubeT0;
+//class MdtTubeFitContainer;
+
+
+
+/**@class T0CalibrationMT
+ Implementation of a T0 calibration using the MT approach.
+ Copied from T0CalibrationClassic
+ @author Felix.Rauscher@Physik.Uni-Muenchen.De
+ @date June 2005
+*/
+class T0CalibrationMT : public IMdtCalibration {
+
+ public:
+
+ /** constructor
+ * @param[in] name of the region/chamber to be calibrated
+ * @param[in] pointer to settings vector
+ * @param[in] sorting criteria (TUBE, CHAMBER, MULTILAYER...) default is by TUBE
+ */
+
+ T0CalibrationMT( std::string name, const T0MTSettings * settings, const std::vector<int> &sort_by, const std::vector<int> &adc_sort_by);
+
+ /** destructor */
+ ~T0CalibrationMT();
+
+ bool handleSegment( MuonCalibSegment& seg ); //!< fill tube spectra
+
+ void setInput( const IMdtCalibrationOutput* input ); //!< unused
+
+ bool analyse(); //!< extract parameters from spectra
+
+ bool converged() const; //!< return m_converged
+
+ /** @return the calibration results
+ */
+ const IMdtCalibrationOutput* getResults() const;
+
+const IMdtCalibrationOutput* analyseSegments( const std::vector<MuonCalibSegment*> & segs );//!< new interface function
+
+ private:
+
+
+ T0MTHistos * getHistos(const MuonFixedId & idtube, unsigned int nr); //!< retrieve pointer for tube idtube histograms
+ ADCMTHistos * getADCHistos(const MuonFixedId & idtube, unsigned int nr); //!< retrieve pointer for tube idtube histograms
+bool analyse_tdc(const int & nr, std::map<int, MdtTubeFitContainer::SingleTubeFit> & full, std::map<int, MdtTubeFitContainer::SingleTubeCalib> & st, std::map<int, std::string> & fit_by_map);
+bool analyse_adc(const int & nr, std::map<int, MdtTubeFitContainer::SingleTubeFit> & full, std::map<int, MdtTubeFitContainer::SingleTubeCalib> & st);
+
+ void doTimeFit(T0MTHistos * T0h, const std::set<MuonFixedId> & tube_ids, std::map<int, MdtTubeFitContainer::SingleTubeFit> & fim, std::map<int, MdtTubeFitContainer::SingleTubeCalib> & stcm, std::map<int, std::string> & fit_by_map, const std::string & fit_by); //!< fit time spectrum
+
+ void doAdcFit(ADCMTHistos * T0h, const std::set<MuonFixedId> & tube_ids, std::map<int, MdtTubeFitContainer::SingleTubeFit> & fim, std::map<int, MdtTubeFitContainer::SingleTubeCalib> & stcm); //!< fit adc spectrum
+
+
+ const T0MTSettings* m_settings; //!< pointer to the settings
+ bool m_converged; //!< convergence status
+ std::string m_name; //!< calibration region name
+ int m_currentItnum; //!< current iteration (always 1?)
+ TFile * p_file; //!< pointer to the histogram file
+ TDirectory * m_regiondir; //!< pointer to the ROOT directory
+
+ std::vector<std::map<HistogramId, T0MTHistos*> >p_histos; //!< vector of pointers tube histograms
+ std::vector<std::map<HistogramId, ADCMTHistos*> >p_adc_histos; //!< vector of pointers tube histograms
+ std::vector<std::map<HistogramId, std::set<MuonFixedId> > > m_tube_ids;
+ std::vector<std::map<HistogramId, std::set<MuonFixedId> > > m_adc_tube_ids;
+ std::map<int, int> m_nhits_per_tube; //!<number of hits per tube
+ std::map<NtupleStationId, MdtTubeFitContainer *> m_result; //!< tube constants
+ std::map<NtupleStationId, MdtRelativeTubeT0> m_rel_tube_t0s;
+ const std::vector<int> &m_sort_by;
+ const std::vector<int> &m_adc_sort_by;
+ bool m_delete_settings;
+};
+}
+
+#endif
+
diff --git a/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibT0/MdtCalibT0/T0CalibrationOutput.h b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibT0/MdtCalibT0/T0CalibrationOutput.h
new file mode 100644
index 0000000000000000000000000000000000000000..48e820aee1c2fc6b9586cf47d97a071469cdb80a
--- /dev/null
+++ b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibT0/MdtCalibT0/T0CalibrationOutput.h
@@ -0,0 +1,44 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#ifndef T0CALIBRATIONOUTPUT_H
+#define T0CALIBRATIONOUTPUT_H
+
+#include "MdtCalibInterfaces/IMdtCalibrationOutput.h"
+#include "MuonCalibStandAloneBase/NtupleStationId.h"
+#include <map>
+
+namespace MuonCalib{
+class MdtTubeFitContainer;
+/**
+@class T0CalibrationOutput
+class for the communication of the results of T0 calibration algorithms
+*/
+class T0CalibrationOutput : public IMdtCalibrationOutput {
+ public:
+ /** constructor (from a pointer to an MdtTubeFitContainer) */
+ T0CalibrationOutput( MdtTubeFitContainer* t0Vec ) :
+ IMdtCalibrationOutput("T0CalibrationOutput"), m_tubeConstants(t0Vec) {}
+ T0CalibrationOutput(const std::map<NtupleStationId, MdtTubeFitContainer*> & t0Map ) :
+ IMdtCalibrationOutput("T0CalibrationOutput"), m_tubeConstants(NULL), m_tubeConstants_map(t0Map) {}
+
+ /** @return the pointer to the MdtTubeFitContainer with the fit results */
+ MdtTubeFitContainer* t0s() const { return m_tubeConstants; }
+ std::map<NtupleStationId, MdtTubeFitContainer*> & GetMap()
+ {
+ return m_tubeConstants_map;
+ }
+ const std::map<NtupleStationId, MdtTubeFitContainer*> & GetMap() const
+ {
+ return m_tubeConstants_map;
+ }
+
+ private:
+ // pointer to a MdtTubeFitContainer instance
+ MdtTubeFitContainer* m_tubeConstants;
+ std::map<NtupleStationId, MdtTubeFitContainer*> m_tubeConstants_map;
+
+};
+}
+#endif
diff --git a/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibT0/MdtCalibT0/T0MTHistos.h b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibT0/MdtCalibT0/T0MTHistos.h
new file mode 100644
index 0000000000000000000000000000000000000000..9265f77769b3f596cc89c55d00e0e7bc7a3f28ce
--- /dev/null
+++ b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibT0/MdtCalibT0/T0MTHistos.h
@@ -0,0 +1,199 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+// T0MTHistos.h
+// Header file for class T0MTHistos
+///////////////////////////////////////////////////////////////////
+// (c) ATLAS Detector software
+///////////////////////////////////////////////////////////////////
+//
+///////////////////////////////////////////////////////////////////
+
+#ifndef T0MTHISTOS_H
+#define T0MTHISTOS_H
+
+//c - c++
+#include <iostream>
+
+//root
+#include "TH1.h"
+//#include "TF1.h"
+//#include "TDirectory.h"
+//class TH1F;
+class TF1;
+class TDirectory;
+
+//this
+#include "T0MTSettings.h"
+
+
+namespace MuonCalib {
+
+/**@class T0MTHistos
+ Histogram and fitter class for drift time and pulsehight spectra
+ The rising slope is fitted by a fermi-function: @f$f(t)=r_{u,0} + \frac{A}{1+exp((t_0 - t)/T_0}@f$.
+ The falling slope is fitted by @f$g(t)=r_{u,max} + \frac{e(t,a,b)}{1+exp((t_{max} - t)/T_{max}}@f$ where @f$e(t, a, b) = a e^{bt}@f$.
+ A pattern recognition determines the fit ranges and the parameters @f$r_u@f$, @f$A@f$, @f$a@f$ and @f$b@f$.
+ @author Felix.Rauscher@Physik.Uni-Muenchen.De
+ @date June 2006
+
+}
+
+*/
+class T0MTHistos
+ {
+ public:
+//---------------------------constructor----------------------------------------
+ /** Default Constructor */
+ inline T0MTHistos() : m_time(NULL),
+ m_id(-1),
+ m_t0_fermi(NULL),
+ m_t0_ok(false),
+ m_status_code(1),
+ m_tmax_fermi(NULL),
+ m_tmax_ok(false),
+ dir(NULL)
+ {
+ m_chi2=9e9;
+ }
+ /** Initializing constructor
+ @param id tube id
+ @param settings t0-fit settings: settings will be asked about histogram binning
+
+ */
+ inline T0MTHistos(int id, const T0MTSettings * settings, const char *hname=NULL):
+ m_t0_fermi(NULL),
+ m_tmax_fermi(NULL)
+ {
+ Initialize(id, settings, hname);
+ m_status_code=99;
+ m_chi2=9e9;
+ }
+//---------------------------static constants-----------------------------------
+ //! number of parameters in t0 fit
+ static const int N_T0_FIT_PAR=4;
+ //! parameter numbers in t0 fit
+ static const int T0_PAR_NR_T0=0,
+ T0_PAR_NR_T=1,
+ T0_PAR_NR_BACK=2,
+ T0_PAR_NR_A=3;
+ //! number of parameters for tmax fit
+ static const int N_TMAX_FIT_PAR=6;
+ //! parameters numbers for tmax fit
+ static const int TMAX_PAR_NR_TMAX=0,
+ TMAX_PAR_NR_T=1,
+ TMAX_PAR_NR_BACK=2,
+ TMAX_PAR_NR_A=3,
+ TMAX_PAR_NR_B=4,
+ TMAX_PAR_NR_T0=5;
+//---------------------------public member functions----------------------------
+ /** Initialize class
+ @param id tube id
+ @param settings t0-fit settings: settings will be asked about histogram binning
+
+ */
+ void Initialize(int id, const T0MTSettings * settings, const char * hname=NULL);
+
+ /** get drift time spectrum */
+ inline TH1F * GetTSpec() {return m_time;}
+
+
+ /** set the pointer of the drift-time spectrum to an existing spectrum.
+ This is for testapps
+ @param id tube id
+ @param spec Pointer to an existing spectrum
+
+ */
+ void SetTSpec(int id, TH1F *spec, const T0MTSettings * settings, bool copy_spec=true);
+
+ /** fill drift time spectrum
+ @param t drift time
+
+ */
+ void FillT(double t);
+
+
+ /** Perform t0-fit
+ Returns true if fit is successfull
+
+ */
+ bool FitT0();
+
+ /** Performs tmax-fit
+ Returns true if fit is successfull
+
+ */
+ bool FitTmax();
+
+ /** return tube id
+
+ */
+ inline int Id() const {return m_id;}
+
+ /** returns true if t0-fit was successfull */
+ inline bool T0Ok() const {return m_t0_ok;}
+ /** returns status code - the status code applies only to the t0 fit*/
+ inline int StatusCode() const {return m_status_code;}
+ /** returns function fitted to the riding edge of the spectrum */
+ inline const TF1 * GetT0Function() const {return m_t0_fermi;}
+ /** returns function fitted to the riding edge of the spectrum */
+ inline TF1 * GetT0Function() {return m_t0_fermi;}
+
+ /** returns true if tmax-fir was successfull */
+ inline bool TmaxOk() const {return m_tmax_ok;}
+
+ /** returns function fitted to the riding edge of the spectrum */
+ inline const TF1 * GetTMaxFunction() const {return m_tmax_fermi;}
+ /** returns function fitted to the riding edge of the spectrum */
+ inline TF1 * GetTMaxFunctionNC() {return m_tmax_fermi;}
+ /** returns t0 chi2*/
+ inline const double & T0Chi2() const
+ {
+ return m_chi2;
+ }
+
+ private:
+//---------------------------private data members-------------------------------
+ //! time spectrum
+ TH1F *m_time;
+ //! tube id;
+ int m_id;
+ //! function fitted to the riding edghe of the spectrum
+ TF1 *m_t0_fermi;
+ //! is true if t0 fit was successful
+ bool m_t0_ok;
+ //! status code for t0 fit (0 ok, 1 not fitted, 2 low statistics, 3 failed)
+ int m_status_code;
+ //! function fitted to the falling edge of the spectrum
+ TF1 *m_tmax_fermi;
+ //! is true if tmax fit was successful
+ double m_tmax_ok;
+ //! TDirectory where debug and result histograms are stored
+ TDirectory *dir;
+ //! Pointer to settings class
+ const T0MTSettings *m_settings;
+ //! chi2/NDF value
+ double m_chi2;
+ /** normal t0 fit */
+ bool NormalFit();
+ /** try to get better start values from a scrambled histogram*/
+ bool T0Scramble();
+ /** top chi2 calculation */
+ void TopChi2();
+ /** top slicing metyhod */
+ void TopSlicing();
+ class Slice
+ {
+ public:
+ double chi_2;
+ int min_bin;
+ int max_bin;
+ int n_bins;
+ };
+
+ };
+
+}
+
+#endif
diff --git a/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibT0/MdtCalibT0/T0MTSettings.h b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibT0/MdtCalibT0/T0MTSettings.h
new file mode 100644
index 0000000000000000000000000000000000000000..c202e9163a0813f810581ffd757a4db4d43aa33a
--- /dev/null
+++ b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibT0/MdtCalibT0/T0MTSettings.h
@@ -0,0 +1,118 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+// T0MTSettings.h
+// Header file for class T0MTSettings
+///////////////////////////////////////////////////////////////////
+// (c) ATLAS Detector software
+///////////////////////////////////////////////////////////////////
+//
+///////////////////////////////////////////////////////////////////
+
+#ifndef T0MTSETTINGS_H
+#define T0MTSETTINGS_H
+
+#include "T0MTSettingsT0.h"
+#include "T0MTSettingsTMax.h"
+
+namespace MuonCalib {
+
+
+
+/**@class T0MTSettings
+ Settings for the T0 calibration (histogram booking and fitting parameters)
+ Parameters for pattern recognition
+ @author Felix.Rauscher@Physik.Uni-Muenchen.De
+ @date Mar 2006
+
+*/
+class T0MTSettings{
+ public:
+ //===========================constructor=======================================
+ /** default constructor*/
+ T0MTSettings();
+//============================public emmber functions===========================
+ /** @name Histogram binning
+ */
+ //@{
+ /** Number of bins for ADC histogram and range*/
+ inline const int & NBinsADC() const {return m_n_bins_adc;}
+ inline int & NBinsADC() {return m_n_bins_adc;}
+ inline const double & ADCMin() const {return m_adc_min;}
+ inline double & ADCMin() {return m_adc_min;}
+ inline const double & ADCMax() const {return m_adc_max;}
+ inline double & ADCMax() {return m_adc_max;}
+ /** Number of bins for time histogram and range*/
+ inline const int & NBinsTime() const {return m_n_bins_time;}
+ inline int & NBinsTime() {return m_n_bins_time;}
+ inline const double & TimeMin() const {return m_time_min;}
+ inline double & TimeMin() {return m_time_min;}
+ inline const double & TimeMax() const {return m_time_max;}
+ inline double & TimeMax() {return m_time_max;}
+ //@}
+ /** @name Debug settings
+ */
+ //@{
+ /** If set to true for every tube a TDirectory will be created. It will contain Graph created from the VariableBinwidthHistograms usded in the Pattern recognition algorithms
+ */
+ inline const bool & DrawDebugGraphs() const {return m_draw_debug_graphs;}
+ inline bool & DrawDebugGraphs() {return m_draw_debug_graphs;}
+ /** If set to true the fitted functions are added to the histograms */
+ inline const bool &AddFitfun() const {return m_add_fitfun;}
+ inline bool &AddFitfun() {return m_add_fitfun;}
+ inline bool &AddFitfun(bool add_fitfun) {m_add_fitfun = add_fitfun;return m_add_fitfun;}
+ /** verbose level
+ 0: no output
+ 1: Fitter output
+ */
+ inline const int &VerboseLevel() const {return m_verbose_level;}
+ inline int &VerboseLevel() {return m_verbose_level;}
+
+ //@}
+ /** get settings for the t0-fit */
+ inline const T0MTSettingsT0 * T0Settings() const {return &m_t0_settings;}
+ inline T0MTSettingsT0 * T0Settings() {return &m_t0_settings;}
+ /** get settings for the tmax-fit*/
+ inline const T0MTSettingsTMax * TMaxSettings() const {return &m_tmax_settings;}
+ inline T0MTSettingsTMax * TMaxSettings() {return &m_tmax_settings;}
+ /** if true fit drift time spectrum*/
+ inline const bool & FitTime() const {return m_fit_time;}
+ inline bool & FitTime() {return m_fit_time;}
+ /** if true fit adc spectrum*/
+ inline const bool & FitADC() const {return m_fit_time;}
+ inline bool & FitADC() {return m_fit_time;}
+ /* minimum number of entries per fit */
+ inline const int & MinEntriesTime() const {return m_min_entries_time;}
+ inline int & MinEntriesTime() {return m_min_entries_time;}
+ inline const int & MinEntriesADC() const {return m_min_entries_adc;}
+ inline int & MinEntriesADC() {return m_min_entries_adc;}
+
+
+private:
+ //===========================private data members=============================
+//----------------------------histogram binning--------------------------------
+ //! ADC bhistogram binning
+ int m_n_bins_adc;
+ double m_adc_min, m_adc_max;
+ //! drift time binning
+ int m_n_bins_time;
+ double m_time_min, m_time_max;
+//---------------------------debug settings-------------------------------------
+ //! if set to true for each debug graphs for the pattern recognition are drawn
+ bool m_draw_debug_graphs;
+ //! if set to true fitted functions are added to hostograms
+ bool m_add_fitfun;
+ //! verbose level
+ int m_verbose_level;
+//---------------------------t0 - tmax - settings-------------------------------
+ bool m_fit_time;
+ bool m_fit_adc;
+ int m_min_entries_time, m_min_entries_adc;
+ T0MTSettingsT0 m_t0_settings;
+ T0MTSettingsTMax m_tmax_settings;
+};
+
+}
+
+#endif
diff --git a/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibT0/MdtCalibT0/T0MTSettingsT0.h b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibT0/MdtCalibT0/T0MTSettingsT0.h
new file mode 100644
index 0000000000000000000000000000000000000000..04092fce540a048dca18996827a9f5418059c78b
--- /dev/null
+++ b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibT0/MdtCalibT0/T0MTSettingsT0.h
@@ -0,0 +1,95 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+// MTT0SettingsT0.h
+// Header file for class MTSettings
+///////////////////////////////////////////////////////////////////
+// (c) ATLAS Detector software
+///////////////////////////////////////////////////////////////////
+//
+///////////////////////////////////////////////////////////////////
+
+#ifndef MTT0SETTINGST0_H
+#define MTT0SETTINGST0_H
+
+/**@class T0MTSettingsT0
+ Settings for the t0 fit to the drift time spectrum
+ @author Felix.Rauscher@Physik.Uni-Muenchen.De
+ @date Mar 2006
+*/
+
+namespace MuonCalib {
+
+class T0MTSettingsT0
+ {
+ public:
+//===========================constructor========================================
+ /** default constructor*/
+ T0MTSettingsT0();
+//===========================public member functions============================
+ /** @name Parameters for pattern recognition
+ Note: These parameters default to well proven parameters. The algorithm should work an all spectra with these. The only parameter which is likely to be adapted is MinBackgroundBins()
+ */
+ //@{
+ /** Number of hits per histogram bin for the VariableBinwidthHistogram
+ The number is given relative to the maximum bin content of the time-Spectrum
+ This number must be >1
+ */
+ inline const double & VBHBinContent() const {return m_vbh_bin_content_rel;}
+ inline double &VBHBinContent() {return m_vbh_bin_content_rel;}
+ /** Maximum bin width for the VariableBinwidthHistogram
+ The bins will not be wider than this even if it means that they will be underfull */
+ inline const double & MaxBinWidth() const {return m_max_bin_width;}
+
+ inline double & MaxBinWidth() {return m_max_bin_width;}
+ /** The minimum width of the region for the background estimation.
+ This width is given in number of bins of the time spectrum.
+ If the pattern recognition selects a smaller region, the lower boarder of the region will be moved to the first bin of the time spectrum.
+
+ NOTE: This should be the approximate width of the trigger matching window before . If this parmeter is to big, the background may be underestimated (with verry little effect on ).
+ */
+ inline const int & MinBackgroundBins() const {return m_min_background_bins;}
+ inline int & MinBackgroundBins() {return m_min_background_bins;}
+ /** Start value for the T-Parameter in the -fit.*/
+ inline const double & Tstart() const {return m_T_start;}
+ inline double & Tstart() {return m_T_start;}
+ /** If true use only the top part of the function for the chi2 calculation*/
+ inline const bool & UseTopChi2() const {return m_use_top_chi2;}
+ inline bool & UseTopChi2(){return m_use_top_chi2;}
+ /** the chi2 threshold at which the scrambling method is used */
+ inline const double & ScrambleThreshold() const {return m_scramble_threshold;}
+ inline double & ScrambleThreshold() {return m_scramble_threshold;}
+ inline double & ScrambleThreshold(const double & scramble_threshold) {m_scramble_threshold = scramble_threshold; return m_scramble_threshold;}
+ /** the chi2 threshold at which the slicing method is used */
+ inline const double & SlicingThreshold() const {return m_slicing_chi2;}
+ inline double & SlicingThreshold() {return m_slicing_chi2;}
+ inline double & SlicingThreshold(const double & slicing_threshold) {m_slicing_chi2 = slicing_threshold;return m_slicing_chi2;}
+ //@}
+ inline const bool & CorrectRelT0s() const { return m_correct_rel_t0s;}
+ inline bool & CorrectRelT0s() { return m_correct_rel_t0s;}
+
+ private:
+ //===========================private data members=============================
+ //! relative bin content for t0-VBH
+ double m_vbh_bin_content_rel;
+ //! maximum bin width for -VBH
+ double m_max_bin_width;
+ //! minimum number of bins for -background estimation
+ int m_min_background_bins;
+ //! T start value
+ double m_T_start;
+ //! use top chi2 instead of total chi2
+ bool m_use_top_chi2;
+ //! chi2 limit for scrambling
+ double m_scramble_threshold;
+ //! chi2 limit for slicing method
+ double m_slicing_chi2;
+ //! correct for relative t0s
+ bool m_correct_rel_t0s;
+
+ };
+
+}
+
+#endif
diff --git a/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibT0/MdtCalibT0/T0MTSettingsTMax.h b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibT0/MdtCalibT0/T0MTSettingsTMax.h
new file mode 100644
index 0000000000000000000000000000000000000000..cff7425fa92d9bf0c49f068638d23fab01d86e7b
--- /dev/null
+++ b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibT0/MdtCalibT0/T0MTSettingsTMax.h
@@ -0,0 +1,100 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+// MTTMaxSettingsTMax.h
+// Header file for class MTSettings
+///////////////////////////////////////////////////////////////////
+// (c) ATLAS Detector software
+///////////////////////////////////////////////////////////////////
+//
+///////////////////////////////////////////////////////////////////
+
+#ifndef MTTMAXSETTINGSTMAX_H
+#define MTTMAXSETTINGSTMAX_H
+
+/**@class T0MTSettingsTMax
+ Settings for the TMax fit to the drift time spectrum
+ @author Felix.Rauscher@Physik.Uni-Muenchen.De
+ @date Mar 2006
+*/
+
+namespace MuonCalib {
+
+class T0MTSettingsTMax
+ {
+ public:
+//===========================constructor========================================
+ /** default constructor*/
+ T0MTSettingsTMax();
+//===========================public member functions============================
+ /** @name Parameters for tmax-pattern recognition
+ Note: These parameters default to well proven parameters. The algorithm should work an all spectra with these. The only parameter which is likely to be adapted is TMaxMinBackgroundBins()
+ */
+ //@{
+ /** Number of hits per histogram bin for the VariableBinwidthHistogram
+ The number is given relative to the maximum bin content of the time-Spectrum
+ This number must be >1
+ */
+ inline const double & VBHBinContent() const {return m_vbh_bin_content_rel;}
+ inline double &VBHBinContent() {return m_vbh_bin_content_rel;}
+ /** Maximum bin width for the VariableBinwidthHistogram
+ The bins will not be wider than this even if it means that they will be underfull */
+ inline const double & MaxBinwidth() const {return m_max_bin_width;}
+
+ inline double & MaxBinwidth() {return m_max_bin_width;}
+ /** For the pattern-recognition only a part of the drift time spectrum is used for creating the VariableBinwidthHistogram. This parameter gives the lower limit for this region minus t0 */
+ inline const double & VBHLow() const {return m_vbh_low;}
+ inline double & VBHLow() {return m_vbh_low;}
+ /** The range of the falling edge of the spectrum is searched in a region after t0. This range is given here. The range defaults to a verry large area which should be ok for all spectra
+ */
+ inline const double & EndMin() const {return m_end_min;}
+ inline double & EndMin() {return m_end_min;}
+ inline const double & EndMax() const {return m_end_max;}
+ inline double & EndMax() {return m_end_max;}
+ /** distance between the detected falling edge and the begin of the background region
+ */
+ inline const double & DistBackground() const {return m_dist_background;}
+ inline double & DistBackground() {return m_dist_background;}
+ /** The minimum width of the region for the background estimation.
+ This width is given in number of bins of the time spectrum.
+ If the pattern recognition selects a smaller region, the upper boarder of the region will be moved to the last bin of the time spectrum.
+
+ NOTE: This should be the approximate width of the trigger matching window after . If this parmeter is to big, the background may be underestimated (with verry little effect on t_max).
+ */
+ inline const int & MinBackgroundBins() const {return m_min_background_bins;}
+ inline int & MinBackgroundBins() {return m_min_background_bins;}
+ /** Distance of the a/b region from the detected falling edge
+ */
+ inline const double & DistAB() const {return m_dist_ab;}
+ inline double & DistAB() {return m_dist_ab;}
+ /** Width of the region in which the parameters a and b are estimated
+ */
+ inline const double & WidthAB() const {return m_width_ab;}
+ inline double & WidthAB() {return m_width_ab;}
+ //@}
+private:
+ //===========================private data members=============================
+ //! relative bin content for tmax-VBH
+ double m_vbh_bin_content_rel;
+ //! maximum bin width for tmax-VBH
+ double m_max_bin_width;
+ //! lower edge of tmax-VBH relative to t0
+ double m_vbh_low;
+ //! range in which the falling edge is searched
+ double m_end_min, m_end_max;
+ //! distance between detected falling edge and begin of background region
+ double m_dist_background;
+ //! minimum number of bins for background region
+ int m_min_background_bins;
+ //! distance of a/b region from detected falling edge
+ double m_dist_ab;
+ //! width of a/b region
+ double m_width_ab;
+ //! start value for T
+ double m_start_t;
+ };
+
+}
+
+#endif
diff --git a/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibT0/MdtCalibT0/VariableBinwidthHistogram.h b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibT0/MdtCalibT0/VariableBinwidthHistogram.h
new file mode 100644
index 0000000000000000000000000000000000000000..a62657fa6eb385b4884745f1b5489de1e92634a4
--- /dev/null
+++ b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibT0/MdtCalibT0/VariableBinwidthHistogram.h
@@ -0,0 +1,157 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+// VariableBinwidthHistogram.h
+// Header file for class VariableBinwidthHistogram
+///////////////////////////////////////////////////////////////////
+// (c) ATLAS Detector software
+///////////////////////////////////////////////////////////////////
+//
+///////////////////////////////////////////////////////////////////
+
+
+#ifndef VARIABLEBINWIDTHHISTOGRAM_H
+#define VARIABLEBINWIDTHHISTOGRAM_H
+
+// c - c++
+#include <algorithm>
+#include <iostream>
+#include <vector>
+
+// ROOT
+//#include "TH1.h"
+//#include "TGraph.h"
+
+class TH1F;
+class TGraph;
+
+//this
+#include "VariableBinwidthHistogramBin.h"
+
+
+namespace MuonCalib {
+
+/** @class VariableBinwidthHistogram
+ A histogram where every bin has the same number of entries. The density is
+ represented by the bin width.
+ @author Felix.Rauscher@Physik.Uni-Muenchen.De
+ @date February 2006
+
+*/
+class VariableBinwidthHistogram {
+ public:
+//----------------------------constructors--------------------------------------
+/** Default constructor */
+ inline VariableBinwidthHistogram() : m_binc(0), m_max_bin_width(0), m_error(false), m_sorted(false) {}
+
+/** destructor */
+ inline ~VariableBinwidthHistogram()
+ {
+ for(unsigned int i=0; i<m_bins.size(); i++)
+ {
+ delete m_bins[i];
+ }
+ }
+
+//-----------------------------public member functinos--------------------------
+
+/** Initialize with new input histogram
+ Returns on error false
+ @param hist the input histogram
+ @param binc_r number of entries per bin / max number of entries in hist
+ @param max_bin_width The binwidth will not exceed the maximum_bin_width
+ @param min_x only the range between min_x and max_x will be used
+ @param max_x only the range between min_x and max_x will be used
+
+*/
+ bool Initialize(TH1F * hist, double binc_r, double max_bin_width, double min_x=-9e9, double max_x=9e9);
+
+
+
+/** Removes steps that origin in a binning effekt
+ @param width smoothin parameter (typically the binwidth of the input-hisatogram)
+
+*/
+ bool Smooth(double width);
+
+/** Get a bin
+ @param bin The bin index of the bin. bin < VariableBinwidthHistogram :: GetNumberOfBins()
+
+ */
+ inline const VariableBinwidthHistogramBin & GetBin(unsigned int bin) const
+ {
+ return *(m_bins[bin]);
+ }
+
+/** Get a bin sorted by content
+ @param bin The bin index of the bin. bin < VariableBinwidthHistogram :: GetNumberOfBins()
+
+ */
+ inline const VariableBinwidthHistogramBin & GetSortedBin(unsigned int bin)
+ {
+ if(!m_sorted)
+ {
+ sort(m_sort_bins.begin(), m_sort_bins.end());
+ m_sorted=true;
+ }
+ return m_sort_bins[bin].Bin();
+ }
+
+/** Get the number of bins
+
+ */
+ inline unsigned int GetNumberOfBins() const
+ {
+ return m_bins.size();
+ }
+
+/** create density graph - density vs bin center */
+ TGraph *DenistyGraph() const;
+
+/** Plot binwidth graph - binwidth versus bin center */
+ TGraph *BinWidthGraph() const;
+
+/** Plot bin content graph - bin content vs bin center */
+ TGraph *BinContentGraph() const;
+
+/** Plot graph with differential density*/
+ TGraph *DiffDensityGraph() const;
+
+/** Plot graph with differential binwidth*/
+ TGraph *DiffBinwidthGraph() const;
+
+ private:
+//----------------------------private data members------------------------------
+
+ /// vector containing histogram bins
+ std::vector<VariableBinwidthHistogramBin *> m_bins;
+
+ /// bins sorted by content
+ std::vector<VBHBinPtrSrt> m_sort_bins;
+
+ /// number of entries per bin
+ unsigned int m_binc;
+
+ /// maximum bin width
+ double m_max_bin_width;
+
+ ///error flag
+ bool m_error;
+
+ bool m_sorted;
+
+//--------------------------privite member functions----------------------------
+ inline double sign(double val) const
+ {
+ if(val>0.0) return 1.0;
+ if(val<0.0) return -1.0;
+ if(val==0.0) return 0;
+ return 0;
+ }
+ };
+
+}
+
+
+#endif
diff --git a/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibT0/MdtCalibT0/VariableBinwidthHistogramBin.h b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibT0/MdtCalibT0/VariableBinwidthHistogramBin.h
new file mode 100644
index 0000000000000000000000000000000000000000..b3e6d1957c7f1d017797185d244aba04d4cc3489
--- /dev/null
+++ b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibT0/MdtCalibT0/VariableBinwidthHistogramBin.h
@@ -0,0 +1,194 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+// VariableBinwidthHistogramBin.h
+// Header file for class VariableBinwidthHistogram
+///////////////////////////////////////////////////////////////////
+// (c) ATLAS Detector software
+///////////////////////////////////////////////////////////////////
+//
+///////////////////////////////////////////////////////////////////
+
+
+#ifndef VARIABLEBINWIDTHHISTOGRAMBIN_H
+#define VARIABLEBINWIDTHHISTOGRAMBIN_H
+
+// c - c++
+#include <cstdlib>
+#include <iostream>
+#include <vector>
+
+namespace MuonCalib {
+
+/** @class VariableBinwidthHistogramBin
+ A bin of a VariableBinwidthHistogram
+ @author Felix.Rauscher@Physik.Uni-Muenchen.De
+ @date February 2006
+
+*/
+class VariableBinwidthHistogramBin {
+ public:
+//----------------------------------constructor---------------------------------
+
+/** Default constructor*/
+ inline VariableBinwidthHistogramBin(): m_center(0.0), m_width(0.0), m_content(0) {}
+
+/** Initializing constructor
+ @param center the center of the bin
+ @param width the width of the bin
+ @param content the number of entries in the bin
+
+*/
+ inline VariableBinwidthHistogramBin(double center, double width, unsigned int content): m_center(center), m_width(width), m_content(content) {}
+
+
+//-----------------------------------public member functions--------------------
+
+/** Set bin
+ @param center the center of the bin
+ @param width the width of the bin
+ @param content the number of entries in the bin
+
+*/
+ inline void Set(double center, double width, unsigned int content)
+ {
+ m_center=center;
+ m_width=width;
+ m_content=content;
+ }
+
+/**set content
+ @param n new bin content
+
+*/
+ inline void SetContent(unsigned int n) {m_content=n;}
+
+/**add to bin content
+ @param n is added
+
+*/
+ inline VariableBinwidthHistogramBin operator += (unsigned int n) {m_content+=n; return *this;}
+/**add to bin content
+ @param n is added
+
+*/
+ inline VariableBinwidthHistogramBin operator + (unsigned int n) {VariableBinwidthHistogramBin ret; ret+=n; return ret;}
+
+/**move right bin boarder
+ @param new_right right bin boarder
+
+*/
+ void MoveRight(double new_right)
+ {
+ double left=Left();
+ if(new_right< left)
+ {
+ std::cerr<<"new right is too small!"<<std::endl;
+ exit(1);
+ }
+ m_width = new_right - left;
+ m_center = 0.5 * (left + new_right);
+ }
+
+/**move left bin boarder
+ @param new_left left bin boarder
+
+*/
+ void MoveLeft(double new_left)
+ {
+ double right=Right();
+ if(new_left>right)
+ {
+ std::cerr<<"new right is too small!"<<std::endl;
+ exit(1);
+ }
+ m_width = right - new_left;
+ m_center = 0.5 * (new_left + right);
+ }
+
+
+/**Get bin center */
+ inline double Center() const {return m_center;}
+
+/**Get width of the bin */
+ inline double Width() const {return m_width;}
+
+/**Get left (lower) bin border */
+ inline double Left() const {return m_center - 0.5 * m_width;}
+
+/**Get right (upper) bin border */
+ inline double Right() const {return m_center + 0.5 * m_width;}
+
+/**Get number of entries in the bin */
+ inline unsigned int Entries() const {return m_content;}
+
+/**Get density=Entries()/Width()*/
+ inline double Density() const {return static_cast<double>(m_content) / m_width;}
+
+
+ private:
+//-----------------------------------------private data members-----------------
+
+ double m_center, m_width;
+ unsigned int m_content;
+ };
+
+/** @class VBHBinPtrSrt
+ A pointer to a VariableBinwidthHistogramBin which supports sorting
+ @author Felix.Rauscher@Physik.Uni-Muenchen.De
+ @date February 2006
+
+*/
+class VBHBinPtrSrt {
+public:
+//-----------------------------------constructor--------------------------------
+/** default constructor */
+ inline VBHBinPtrSrt() {}
+
+/** initializing constructor
+ @param bin Pointer to bin
+
+*/
+ inline VBHBinPtrSrt(VariableBinwidthHistogramBin *bin): m_bin(bin) {}
+
+//--------------------------------------public member functions-----------------
+
+/** initialize
+ @param bin Pointer to bin
+
+*/
+ inline void Initialize(VariableBinwidthHistogramBin *bin) {m_bin=bin;}
+
+/** Get reference to bin
+
+*/
+ inline VariableBinwidthHistogramBin & Bin() {return *m_bin;}
+
+/**Operator > for sorting bins by content
+ @param other Other bin
+
+*/
+ inline bool operator > (const VBHBinPtrSrt &other) const
+ {
+ return m_bin->Width() > other.m_bin->Width();
+ }
+
+/**Operator < for sorting bins by content
+ @param other Other bin
+
+*/
+ inline bool operator < (const VBHBinPtrSrt &other) const
+ {
+ return m_bin->Width() < other.m_bin->Width();
+ }
+private:
+//------------------------------private data members----------------------------
+
+ /// pointer to bin
+ VariableBinwidthHistogramBin *m_bin;
+ };
+
+}
+
+#endif
diff --git a/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibT0/MdtCalibT0/selection.xml b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibT0/MdtCalibT0/selection.xml
new file mode 100644
index 0000000000000000000000000000000000000000..709c823abf91ab8dd11fc151f7655dea8a3dd108
--- /dev/null
+++ b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibT0/MdtCalibT0/selection.xml
@@ -0,0 +1,8 @@
+<lcgdict>
+
+ <class name="MuonCalib::T0MTHistos" />
+ <class name="MuonCalib::T0MTSettings" />
+ <class name="MuonCalib::T0MTSettingsT0" />
+ <class name="MuonCalib::T0MTSettingsTMax" />
+
+</lcgdict>
diff --git a/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibT0/cmt/requirements b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibT0/cmt/requirements
new file mode 100644
index 0000000000000000000000000000000000000000..63aaca313080d0cf578c27fd2cfc7157320dcefa
--- /dev/null
+++ b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibT0/cmt/requirements
@@ -0,0 +1,37 @@
+package MdtCalibT0
+
+author Domizia Orestano
+
+use AtlasPolicy AtlasPolicy-*
+
+private
+apply_tag ROOTMathLibs
+apply_tag ROOTGraphicsLibs
+apply_tag ROOTRooFitLibs
+end_private
+use AtlasROOT AtlasROOT-* External
+
+use MdtCalibInterfaces * MuonSpectrometer/MuonCalib/MdtCalib
+use MdtCalibData * MuonSpectrometer/MuonCalib/MdtCalib
+
+use MuonCalibStandAloneBase * MuonSpectrometer/MuonCalib/MuonCalibStandAlone
+
+#macro MdtCalibT0_linkopts "-lMdtCalibT0 -lcurses -lCore -lCint -lHist -lGraf -lGraf3d -lGpad -lTree -lRint -lPostscript -lMatrix -lPhysics -pthread -lm -ldl -rdynamic "
+
+private
+
+use MuonCalibEventBase * MuonSpectrometer/MuonCalib
+use MuonCalibIdentifier * MuonSpectrometer/MuonCalib
+use MuonCalibStl * MuonSpectrometer/MuonCalib/MuonCalibUtils
+
+
+end_private
+
+library MdtCalibT0 ../src/*.cxx
+apply_pattern installed_library
+
+#macro_append T0Fit_dependencies "MdtCalibT0"
+#application T0Fit ../exe/T0Fit.cxx
+
+use AtlasReflex AtlasReflex-* External
+apply_pattern lcgdict dict=MdtCalibT0 selectionfile=selection.xml headerfiles="../MdtCalibT0/ReflexHeaders.h "
diff --git a/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibT0/doc/mainpage.h b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibT0/doc/mainpage.h
new file mode 100755
index 0000000000000000000000000000000000000000..26ac53fdb550548cd01fb2af0b94fdd148a0677d
--- /dev/null
+++ b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibT0/doc/mainpage.h
@@ -0,0 +1,48 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+/**
+@mainpage MdtCalibT0 Package
+@author Domizia.Orestano@cern.ch
+
+@section MdtCalibT0Intro Introduction
+This package is foreseen to collect the various implementations of Mdt T0 calibration.
+
+@section MdtCalibT0Overview Class Overview
+The MdtCalibT0 package currently contains the following classes common to all the implementations:
+
+- MuonCalib::T0CalibrationOutput: Class for communication between event loop and T0 calibration
+
+The following classes related to "Classic" implementation:
+- MuonCalib::T0ClassicSettings: Parameters needed by the T0 calibrations which should be provided to the T0CalibrationClassic constructor
+- MuonCalib::T0CalibrationClassic: Implementation of a T0 calibration using the classical approach
+- MuonCalib::T0ClassicHistos: Tube histograms used in T0 calibration
+
+The following classes related to "MT" implementation:
+- MuonCalib::T0MTSettings: Parameters needed by the T0 calibrations which should be provided to the T0CalibrationMT constructor. Uses T0MTSettingsT0 and T0MTSettingsTMax.
+- MuonCalib::T0MTSettingsT0: Settings for the T0 fit to the drift time spectrum
+- MuonCalib::T0MTSettingsTMax: Settings for the TMax fit to the drift time spectrum
+- MuonCalib::T0CalibrationMT: Implementation of a T0 calibration using the MT approach
+- MuonCalib::T0MTHistos:Histogram and fitter class for drift time and pulsehight spectra.
+- MuonCalib::MTT0PattternRecognition: Performs the pattern recognition for the MT t0-fit routine. Fir ranges and start values are estimated.
+- MuonCalib::MTTmaxPattternRecognition: Performs the pattern recognition for the MT tmax-fit routine. Fir ranges and start values are estimated.
+- MuonCalib::HistogramId:Identifier class for drift time histograms. The class can be switched to the identification of single tubes, and different types of groupes of tubes.
+- MuonCalib::VariableBinwidthHistogram: a histogram where every bin has the same number of entries. The density is represented by the bin width.
+- MuonCalib::VariableBinwidthHistogramBin: a bin of a VariableBinwidthHistogram.
+
+@ref used_MdtCalibT0
+
+@ref requirements_MdtCalibT0
+*/
+
+/**
+@page used_MdtCalibT0 Used Packages
+@htmlinclude used_packages.html
+*/
+
+/**
+@page requirements_MdtCalibT0 Requirements
+@include requirements
+*/
+
diff --git a/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibT0/exe/T0Fit.cxx b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibT0/exe/T0Fit.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..6e5f2e66bd3b98fadc5727b2862861e06a902d68
--- /dev/null
+++ b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibT0/exe/T0Fit.cxx
@@ -0,0 +1,69 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#include "MdtCalibT0/T0MTHistos.h"
+#include "TROOT.h"
+#include "TFile.h"
+#include "TH1.h"
+#include "TDirectory.h"
+#include "TKey.h"
+#include "TClass.h"
+#include "TObject.h"
+#include "iostream"
+#include "string"
+#include "cstdlib"
+
+using namespace MuonCalib;
+using namespace std;
+
+
+
+int main(int argc, char * argv[])
+ {
+//check command line arguments
+ if(argc!=3)
+ {
+ cerr<<"Usage: "<<argv[0]<<"<input file> <output file>!"<<endl;
+ std::exit(1);
+ }
+
+
+ TROOT wurscht("wurscht", "wurscht");
+ TFile infile(argv[1]);
+ TDirectory *dir=&infile;
+ TFile outfile(argv[2], "RECREATE");
+ T0MTSettings settings;
+ settings.AddFitfun()=true;
+ settings.DrawDebugGraphs()=true;
+ if (dir==NULL)
+ {
+ cerr<<"Cannot find TDirectory 'MT_t0_fitter' in file '"<<argv[1]<<endl;
+ exit(1);
+ }
+//loop on all histograms
+ TIter nextkey(dir->GetListOfKeys());
+ TKey *key;
+ int n=0;
+ while((key=(TKey *) nextkey()))
+ {
+ TObject *obj = key->ReadObj();
+ if (obj->IsA()->InheritsFrom("TH1F"))
+ {
+ TH1F *hist=dynamic_cast<TH1F *>(obj);
+ if (hist==NULL) continue;
+ string hname=hist->GetName();
+ if (hname=="t_spec_Summary") continue;
+// if((int)hname.find("t_spec")<0) continue;
+ hist->GetListOfFunctions()->Clear();
+ T0MTHistos fitter;
+ fitter.SetTSpec(n, hist, &settings);
+ n++;
+ if(fitter.FitT0())
+ {
+ fitter.FitTmax();
+ }
+ }
+ }
+ outfile.Write();
+ }
diff --git a/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibT0/src/ADCMTHistos.cxx b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibT0/src/ADCMTHistos.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..586fc85df9cab368aa1b2d9a0e43681f75dbfbe2
--- /dev/null
+++ b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibT0/src/ADCMTHistos.cxx
@@ -0,0 +1,103 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+using namespace std;
+
+#include "MdtCalibT0/ADCMTHistos.h"
+
+//root
+#include "TH1.h"
+#include "TF1.h"
+#include "TMath.h"
+#include "TDirectory.h"
+#include <cmath>
+
+namespace MuonCalib {
+
+/** skew normal ADC fitting
+*/
+ inline Double_t adcfitf_skewnormal(Double_t *x, Double_t *par){
+ //par[0] = skew gauss norm
+ //par[1] = skew gauss mean (i.e. mu)
+ //par[2] = skew gauss sigma (i.e sigma)
+ //par[3] = skew factor (i.e. alpha)
+ // Numeric constants
+ //Double_t invsq2pi = 0.3989422804014; // (2 pi)^(-1/2)
+ Double_t invsq2pi = 1.0/(TMath::Sqrt(2*3.14159265358979323846264338327950288419716939937510));
+ Double_t twoPi = 6.2831853071795;//2Pi
+
+// Double_t Landau_Value = par[0]*TMath::Landau(x[0],par[1],par[2],kFALSE);
+
+ Double_t delta_value = par[3]/(TMath::Sqrt(1.+par[3]*par[3]));
+ Double_t omega_square = (par[2]*par[2])/(1. - 4.*delta_value*delta_value/(twoPi));
+ Double_t omega_value = TMath::Sqrt(omega_square);
+ Double_t xi_value = par[1] - delta_value*omega_value*2.*invsq2pi;
+ Double_t Gauss_part= (invsq2pi/omega_value)*exp(-((x[0] - xi_value)*(x[0] - xi_value))/(2.0*omega_square));//phi(x)
+
+ Double_t Erf_part=0.5*( 1+TMath::Erf(par[3]*(x[0]-xi_value)/(omega_value)) );
+
+ Double_t SkewNormal_Value = par[0]*2.*Gauss_part*Erf_part;
+
+ Double_t MyGaussFuncValue = SkewNormal_Value;
+
+ return(MyGaussFuncValue);
+
+ }
+
+
+void ADCMTHistos :: FillA(double a)
+ {m_adc->Fill(static_cast<Axis_t>(a));}
+
+
+//////////////////////////////////////////////////////////
+// Initialize(int id, const T0MTSettings & settings) //
+//////////////////////////////////////////////////////////
+
+void ADCMTHistos :: Initialize(int id, const T0MTSettings * settings, const char * hname)
+ {
+ m_settings=settings;
+// if(m_settings->VerboseLevel()>1)
+// {
+ cout<<"ADCMTHistos :: Initialize: called"<<endl;
+// }
+ char buf[100];
+ std::cout<<gDirectory->GetName()<<std::endl;
+ if(hname==NULL)
+ snprintf(buf, 100, "adc_spec_%d", id);
+ else
+ snprintf(buf, 100, "adc_spec_%s", hname);
+ m_adc=new TH1F(buf, "", settings->NBinsADC(), settings->ADCMin(), settings->ADCMax());
+ m_id=id;
+ m_adc_ok=false;
+ }
+
+
+
+bool ADCMTHistos :: FitAdc()
+ {
+/* m_adc_fit = new TF1("adc_fun", "landau(0) + gaus(3)");
+ m_adc->Fit("landau", "Q");
+ TF1 *lan_fit=m_adc->GetFunction("landau");
+ for(int i=0; i<3; i++)
+ m_adc_fit->SetParameter(i, lan_fit->GetParameter(i));
+ m_adc->Fit("adc_fun", "Q");
+*/
+ m_adc_fit = new TF1("adc_fun",adcfitf_skewnormal, 50, 320, 4);
+ Double_t average = m_adc->GetMean();
+ Double_t max=m_adc->GetMaximum();
+// m_adc->SetAxisRange(50,350,"X");
+ m_adc_fit->SetParameters(max,average,40,1); // initialize value
+ m_adc_fit->SetLineColor(kRed);
+ m_adc_fit->SetParNames("Max","Mean", "Sigma","Skew_factor");
+ m_adc_fit->SetParLimits(0,0,1000000);
+ m_adc_fit->SetParLimits(1,100,200);
+ m_adc_fit->SetParLimits(2,0,200);
+ //gStyle->SetOptFit(1011);
+
+ m_adc->Fit("adc_fun","Q+","",50,320);
+ m_adc_ok = true;
+ return true;
+ }
+
+} //namespace MuonCalib
diff --git a/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibT0/src/HistogramId.cxx b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibT0/src/HistogramId.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..28329c4627083abb437afb5dc387b0e9576d9631
--- /dev/null
+++ b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibT0/src/HistogramId.cxx
@@ -0,0 +1,117 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+//c - c++
+#include <algorithm>
+#include "iostream"
+#include "sstream"
+#include "iomanip"
+
+//athena
+#include "MdtCalibT0/HistogramId.h"
+#include "MuonCalibIdentifier/MuonFixedId.h"
+
+namespace MuonCalib {
+
+std::map<std::pair<int, int>, std::string> HistogramId :: s_histogram_names;
+
+
+void HistogramId :: Initialize(const MuonFixedId &id, int sort_by)
+ {
+ m_id.first=sort_by;
+ static const unsigned int kStationNameShift = 24;
+
+ static const unsigned int kStationEtaShift = 19;
+
+ static const unsigned int kStationPhiShift = 13;
+
+
+// Mdt specific code
+
+ static const unsigned int kMdtMultilayerShift = 9;
+
+ static const unsigned int kMdtTubeLayerShift = 7;
+
+// static const unsigned int kMdtTubeShift = 0;
+ switch(sort_by)
+ {
+ case TUBE:
+ {
+ m_id.second = id.getIdInt();
+ break;
+ }
+ case LAYER:
+ {
+ m_id.second = (id.stationNameIndex()<<kStationNameShift) | (id. etaIndex()<<kStationEtaShift) | (id.phiIndex()<<kStationPhiShift) | (id.mdtMultilayerIndex()<<kMdtMultilayerShift) | (id.mdtTubeLayerIndex()<<kMdtTubeLayerShift);
+ break;
+ }
+ case MULTILAYER:
+ {
+ m_id.second = (id.stationNameIndex()<<kStationNameShift) | (id. etaIndex()<<kStationEtaShift) | (id.phiIndex()<<kStationPhiShift) | (id.mdtMultilayerIndex()<<kMdtMultilayerShift);
+// std::cout<<"name"<<id.stationNameIndex()<<std::endl;
+// std::cout<<"eta"<<id.etaIndex()<<std::endl;
+// std::cout<<"ml"<<id.mdtMultilayerIndex()<<std::endl;
+ break;
+ }
+ case CHAMBER:
+ {
+ m_id.second = (id.stationNameIndex()<<kStationNameShift) | (id. etaIndex()<<kStationEtaShift) | (id.phiIndex()<<kStationPhiShift);
+ break;
+ }
+ case MEZZ_CARD:
+ {
+ //first part is like multilayer
+ m_id.second = id.mdtMezzanine();
+ break;
+ }
+ default:
+ {
+ std::cerr<<"HistogramId :: Initialize: sort_by arguemnt is invalid!"<<std::endl;
+ }
+ }
+//create histogram name
+ if(s_histogram_names[m_id] == "")
+ {
+ std::ostringstream os;
+ if(m_id.second==-999999)
+ {
+ os<<"Summary";
+ }
+ else
+ {
+ //chamber name
+ os << id.stationNumberToFixedStationString(id.stationName()) <<"_eta";
+ if(id.eta()<0) os<<"C";
+ else os<<"A";
+ os<<abs(id.eta())<<"_phi"<<id.phi();
+ //multilayer name
+ if(sort_by != CHAMBER)
+ os<<"_ml"<<id.mdtMultilayer();
+ //layer name
+ if(sort_by != CHAMBER && sort_by != MULTILAYER && sort_by!=MEZZ_CARD)
+ os<<"_ly"<<id.mdtTubeLayer();
+ //tube name
+ if(sort_by == TUBE)
+ os<<"_tb"<<id.mdtTube();
+ //mezz-id
+ if(sort_by == MEZZ_CARD)
+ {
+
+ os<<"_mez"<<(id.mdtMezzanine() %100);
+ }
+ //numeric_id
+ os<<"_num"<<m_id.second<<"_"<<sort_by;
+
+ }
+ //store
+ s_histogram_names[m_id] = os.str();
+ std::cout<<sort_by<<" "<<s_histogram_names[m_id]<<std::endl;
+ }
+ }
+
+
+
+
+
+}
diff --git a/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibT0/src/MTT0PatternRecognition.cxx b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibT0/src/MTT0PatternRecognition.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..2b97a9b3cc59333b2d19784bc4e7258bf300dc0a
--- /dev/null
+++ b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibT0/src/MTT0PatternRecognition.cxx
@@ -0,0 +1,167 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+using namespace std;
+
+//this
+#include "MdtCalibT0/MTT0PatternRecognition.h"
+
+//root
+#include "TLine.h"
+#include "TGraph.h"
+#include "TH1.h"
+
+namespace MuonCalib{
+
+
+//////////////////////////////////////////////////////////////////
+// estimate_background(const TH1 &hist, double scale_min) //
+//////////////////////////////////////////////////////////////////
+
+
+bool MTT0PatternRecognition :: estimate_background( TH1F *hist, double scale_min)
+ {
+//Get first bin in input histogram which is not empty. This avoids underestimation of the background rate if the range of the input histogram exceeds the TDC-range
+ int min(-1);
+ for(int i=1; i<=hist->GetNbinsX(); i++)
+ {
+ if(hist->GetBinContent(i)>0)
+ {
+ min=i;
+ break;
+ }
+ }
+ if(min==-1)
+ {
+ cerr<<"MTT0PattternRecognition :: estimate_background: No hits in input histogram!"<<endl;
+ m_error=true;
+ return false;
+ }
+ double maxx=scale_min-40;
+ int max=hist->FindBin(maxx);
+//we want 10 bins minimum for the background estimate
+ if(max-min<m_settings->MinBackgroundBins())
+ {
+ min=max-128;
+ if (min<0)
+ min=0;
+ }
+//if there are still not enough bins for the background estimate we have a problem
+ if(max-min<m_settings->MinBackgroundBins())
+ {
+ cerr<<"MTT0PattternRecognition :: estimate_background: Rising edge is to glose to lower histogram range!"<<endl;
+ m_error=true;
+ return false;
+ }
+//calculate average bin content
+ m_background=0.0;
+ double back_squared=0.0;
+ double n_bins=0.0;
+ double referece_chi2=0.0;
+ std::cout<<"UUUuuuUUU ";
+ for(int i=min; i<max; i++)
+ {
+ n_bins++;
+ m_background+=hist->GetBinContent(i);
+ back_squared+=hist->GetBinContent(i) * hist->GetBinContent(i);
+ if (n_bins==m_settings->MinBackgroundBins())
+ {
+ double bac=m_background/n_bins;
+ referece_chi2 = 2 * (back_squared/n_bins - bac*bac);
+ std::cout<<referece_chi2<<" ";
+ }
+ if (n_bins>m_settings->MinBackgroundBins())
+ {
+ double bac=m_background/n_bins;
+ double chi2=2 * (back_squared/n_bins - bac*bac);
+ std::cout<<chi2<<" ";
+ if (chi2>5*referece_chi2) break;
+ }
+ }
+ std::cout<<std::endl;
+ m_background/=n_bins;
+//store lower edge of fit range
+ m_fit_min=hist->GetBinCenter(min);
+ m_t0_est=0.5*(scale_min + hist->GetBinCenter(max));
+//mark selected range
+ if(m_draw_debug_graph)
+ {
+ (new TLine(hist->GetBinCenter(min), 0, hist->GetBinCenter(min), hist->GetMaximum()))->Write("t0_back_left");
+ (new TLine(hist->GetBinCenter(max), 0, hist->GetBinCenter(max), hist->GetMaximum()))->Write("t0_back_right");
+ }
+
+ return true;
+ }
+
+
+//////////////////////////////////////////
+// estimate_height(const TH1 &hist) //
+//////////////////////////////////////////
+
+double MTT0PatternRecognition :: estimate_height( TH1F *hist)
+ {
+//smooth histogram
+ if(!m_vbh.Smooth(hist->GetBinWidth(1)))
+ {
+ m_error=true;
+ return -1;
+ }
+ if(!m_vbh.Smooth(hist->GetBinWidth(1)/2))
+ {
+ m_error=true;
+ return -1;
+ }
+//draw debug graphs
+ if(m_draw_debug_graph)
+ {
+ m_vbh.DenistyGraph()->Write("t0_density");
+ m_vbh.BinWidthGraph()->Write("t0_binwidth");
+ m_vbh.BinContentGraph()->Write("t0_bin_content");
+ m_vbh.DiffDensityGraph()->Write("t0_diff_density");
+ m_vbh.DiffDensityGraph()->Write("t0_diffbinwidth");
+ }
+//get the range in which the smallest, and second smallest bins are
+ double lower=m_vbh.GetSortedBin(0).Width();
+ double left(9e9), right(-9e9);
+// bool smallest_found(false);
+ for(unsigned int i=0; i<m_vbh.GetNumberOfBins(); i++)
+ {
+ if(m_vbh.GetSortedBin(i).Left()<left) left = m_vbh.GetSortedBin(i).Left();
+ if(m_vbh.GetSortedBin(i).Right()>right) right = m_vbh.GetSortedBin(i).Right();
+//NOTE: For greater numerical stability m_vbh.GetSortedBin(i).Width() is considered greater only if it is greater by hist->GetBinWidth(1) / 100.0. I have seen differences in Binwidth in the order of 1e-13. On the other hand, due to the binning of the input histogram, e real difference in the binwidth of the VBH n*hist->GetBinWidth(1).
+ if((m_vbh.GetSortedBin(i).Width() - lower) > hist->GetBinWidth(1) / 100.0)
+ {
+ //reqire minimum distance between bins, and minimum number if significant bins
+ if(right-left>50 && i>20) break;
+ lower=m_vbh.GetSortedBin(i).Width();
+ }
+ }
+//mark selected range
+ if(m_draw_debug_graph)
+ {
+ (new TLine(left, 0, left, hist->GetMaximum()))->Write("t0_height_left");
+ (new TLine(right, 0, right, hist->GetMaximum()))->Write("t0_height_right");
+ }
+//calcularte mean bin content of input histogram in range
+ int bleft(hist->FindBin(left)), bright(hist->FindBin(right));
+ double nbins=0.0;
+ m_height=0.0;
+ for(int i=bleft; i<=bright; i++)
+ {
+ m_height+=hist->GetBinContent(i);
+ nbins++;
+ }
+ if(nbins<1)
+ {
+ cerr<<"top region is too small!"<<endl;
+ cout<<"left="<<left<<" right="<<right<<endl;
+ m_error=true;
+ return -1;
+ }
+ m_height/=nbins;
+ m_fit_max=right;
+ return left;
+ }
+
+}
diff --git a/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibT0/src/MTTmaxPatternRecognition.cxx b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibT0/src/MTTmaxPatternRecognition.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..efc6ee6c4c75435828a405fe9c124610ca4f65b0
--- /dev/null
+++ b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibT0/src/MTTmaxPatternRecognition.cxx
@@ -0,0 +1,140 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+using namespace std;
+
+#include "MdtCalibT0/MTTmaxPatternRecognition.h"
+
+//root
+#include "TF1.h"
+#include "TLine.h"
+#include "TGraph.h"
+#include "TH1.h"
+#include <cmath>
+
+namespace MuonCalib{
+
+
+
+inline Double_t myexp(Double_t *x, Double_t *p)
+ {
+ return exp(p[0] + p[1] * (x[0] - p[3])) + p[2];
+ }
+
+bool MTTmaxPatternRecognition :: Initialize( TH1F *hist, double t0, const T0MTSettings * settings)
+ {
+ m_settings=settings->TMaxSettings();
+ m_draw_debug_graph=settings->DrawDebugGraphs();
+ m_error=false;
+ m_t0=t0;
+ m_error=!m_vbh.Initialize(hist, m_settings->VBHBinContent(), m_settings->MaxBinwidth(), t0 + m_settings->VBHLow());
+ if(m_error || m_vbh.GetNumberOfBins()<2)
+ {
+ cerr<<"MTTmaxPatternRecognition :: Initialize: INitialization of VBH failed!"<<endl;
+ }
+ if(!m_error) estimate_background(hist);
+ if(!m_error) estimate_height(hist);
+ return !m_error;
+ }
+
+
+//////////////////////////////////////////
+// estimate_background( TH1 *hist) //
+//////////////////////////////////////////
+
+bool MTTmaxPatternRecognition :: estimate_background( TH1F *hist)
+ {
+//smooth VBH in several steps
+ for(double i=1.0; i<=4.0; i*=2.0)
+ {
+ for(int j=0; j<2; j++)
+ if(!m_vbh.Smooth(hist->GetBinWidth(1)/i))
+ {
+ m_error=true;
+ return -1.0;
+ }
+ }
+ if(m_draw_debug_graph)
+ {
+ m_vbh.DenistyGraph()->Write("tmax_density");
+ m_vbh.BinWidthGraph()->Write("tmax_binwidth");
+ m_vbh.DiffDensityGraph()->Write("tmax_diffdensity");
+ m_vbh.DiffBinwidthGraph()->Write("tmax_diffbinwidth");
+ }
+//find falling edge - search for steepenst slope in a range of 600 to 800 ns after the t0
+ double peak_falling(0.0), peak_falling_slope(-9e9);
+ for(unsigned int i=0; i<m_vbh.GetNumberOfBins()-1; i++)
+ {
+ const VariableBinwidthHistogramBin &bin1(m_vbh.GetBin(i)), &bin2(m_vbh.GetBin(i+1));
+ //select region
+ if(bin1.Center()<m_t0+m_settings->EndMin() || bin1.Center()>m_t0+m_settings->EndMax()) continue;
+ if(bin2.Width() - bin1.Width() > peak_falling_slope && bin2.Right() < hist->GetBinLowEdge(hist->GetNbinsX()) - 10)
+ {
+ peak_falling=bin1.Right();
+ peak_falling_slope=bin2.Width() - bin1.Width();
+ }
+ }
+// cout<<"Peak falling="<<peak_falling<<endl;
+//check region size
+ int firstbin=hist->FindBin(peak_falling+m_settings->DistBackground());
+ int lastbin=hist->FindBin(m_vbh.GetBin(m_vbh.GetNumberOfBins() - 1).Right());
+ if(lastbin-firstbin<m_settings->MinBackgroundBins()) lastbin=hist->GetNbinsX();
+ if(lastbin-firstbin<m_settings->MinBackgroundBins())
+ {
+ cerr<<"MTTMaxPattternRecognition :: estimate_background: Falling edge is to glose to upper histogram range!"<<endl;
+ m_error=true;
+ return false;
+ }
+//calcul;ate mean
+ m_background=0.0;
+ double n_bins=0.0;
+ if(lastbin>hist->GetNbinsX())lastbin=hist->GetNbinsX();
+ for(int i=firstbin; i<=lastbin; i++)
+ {
+ m_background += hist->GetBinContent(i);
+ n_bins++;
+ }
+//mark selected range
+ if(m_draw_debug_graph)
+ {
+ (new TLine(hist->GetBinCenter(firstbin), 0, hist->GetBinCenter(firstbin), hist->GetMaximum()))->Write("tmax_back_left");
+ (new TLine(hist->GetBinCenter(lastbin), 0, hist->GetBinCenter(lastbin), hist->GetMaximum()))->Write("tmax_back_right");
+ }
+
+ m_background/=n_bins;
+ m_tmax_est=peak_falling;
+ m_fit_max=hist->GetBinCenter(lastbin);
+ return true;
+ }
+
+
+//////////////////////////////////////////////////////////
+// estimate_height( TH1 *hist, double background_min) //
+//////////////////////////////////////////////////////////
+
+bool MTTmaxPatternRecognition :: estimate_height( TH1F *hist)
+ {
+//get start value for fit
+ Double_t left = m_tmax_est - m_settings->DistAB() - m_settings->WidthAB(),
+ right = m_tmax_est - m_settings->DistAB();
+//fit
+ TF1 *fun=new TF1("myexp", myexp, left, right, 4);
+ fun->FixParameter(2, m_background);
+ fun->FixParameter(3, m_t0);
+ hist->Fit("myexp", "Q0+", "", left, right);
+//store parameters
+// TF1 *fun=hist->GetFunction("expo");
+ m_a=fun->GetParameter(0);
+ m_b=fun->GetParameter(1);
+ m_fit_min=m_tmax_est - m_settings->DistAB() - m_settings->WidthAB();
+//mark selected range
+ if(m_draw_debug_graph)
+ {
+ (new TLine(left, 0, left, hist->GetMaximum()))->Write("tmax_height_left");
+ (new TLine(right, 0, right, hist->GetMaximum()))->Write("tmax_height_right");
+ }
+ return true;
+ }
+
+}
diff --git a/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibT0/src/MdtRelativeTubeT0.cxx b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibT0/src/MdtRelativeTubeT0.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..e88e782b1532c6aed91e6f2ecdc0fa7014e986b1
--- /dev/null
+++ b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibT0/src/MdtRelativeTubeT0.cxx
@@ -0,0 +1,101 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#include "MdtRelativeTubeT0.h"
+
+//MuonCalibIdentifier
+#include "MuonCalibIdentifier/MuonFixedId.h"
+
+//MuonCalibEventBase
+#include "MuonCalibEventBase/MdtCalibHitBase.h"
+
+//#include <pair>
+
+namespace MuonCalib {
+
+inline unsigned int get_group_id(const MuonFixedId & id, const MdtRelativeTubeT0::TubeGroup &grp);
+
+void MdtRelativeTubeT0 :: AddHit(const MdtCalibHitBase * hit)
+ {
+ MuonFixedId id(hit->identify());
+ if(m_tube_t0.find(id)==m_tube_t0.end())
+ {
+ if( m_relative_offset.size() ) m_relative_offset.clear();
+ m_tube_t0[id] = hit->tubeT0();
+ }
+ }
+
+double MdtRelativeTubeT0 :: GetRelativeOffset(const MuonFixedId & id, TubeGroup grp) const
+ {
+ if (grp==UNKNOWN) return 0;
+ std::map<TubeGroup, std::map<MuonFixedId, double> > :: const_iterator it=m_relative_offset.find(grp);
+ if(it==m_relative_offset.end())
+ {
+ calculate_relative_t0s(grp);
+ return GetRelativeOffset(id, grp);
+ }
+ std::map<MuonFixedId, double> :: const_iterator it2=it->second.find(id);
+ if(it2==it->second.end())
+ {
+ return 0;
+ }
+ return it2->second;
+ }
+
+inline void MdtRelativeTubeT0 :: calculate_relative_t0s(const TubeGroup & grp) const
+ {
+//calculate mean t0 per group
+ std::map<unsigned int, std::pair<double, int> > mean_t0;
+ for(std::map<MuonFixedId, double> :: const_iterator it=m_tube_t0.begin(); it!=m_tube_t0.end(); it++)
+ {
+ unsigned int grp_id(get_group_id(it->first, grp));
+ std::map<unsigned int, std::pair<double, int> > :: iterator it2=mean_t0.find(grp_id);
+ if(it2==mean_t0.end())
+ {
+ mean_t0[grp_id] = std::pair<double, int>(it->second, 1.);
+ }
+ else
+ {
+ it2->second.first+=it->second;
+ it2->second.second++;
+ }
+ }
+ for(std::map<unsigned int, std::pair<double, int> > :: iterator it=mean_t0.begin(); it!=mean_t0.end(); it++)
+ {
+ it->second.first/=it->second.second;
+ }
+
+//calculate tube offsets
+ std::map<MuonFixedId, double> &offsets(m_relative_offset[grp]);
+ offsets.clear();
+ for (std::map<MuonFixedId, double> :: const_iterator it=m_tube_t0.begin(); it!=m_tube_t0.end(); it++)
+ {
+ unsigned int grp_id(get_group_id(it->first, grp));
+ offsets[it->first] = it->second - mean_t0[grp_id].first;
+ }
+ }
+
+inline unsigned int get_group_id(const MuonFixedId & id, const MdtRelativeTubeT0::TubeGroup &grp)
+ {
+ switch(grp)
+ {
+ case MdtRelativeTubeT0::CHAMBER:
+ return id.mdtChamberId().getIdInt();
+ break;
+ case MdtRelativeTubeT0::MULTILAYER:
+ return id.mdtMultilayerId().getIdInt();
+ break;
+ case MdtRelativeTubeT0::LAYER:
+ return 4*(id.mdtMultilayer() -1 ) + (id.mdtTubeLayer() -1);
+ break;
+ case MdtRelativeTubeT0::MEZZ_CARD:
+ return id.mdtMezzanine();
+ break;
+ case MdtRelativeTubeT0::UNKNOWN:
+ return 0;
+ }
+ return 0;
+ }
+
+} //namespace MuonCalib
diff --git a/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibT0/src/MdtRelativeTubeT0.h b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibT0/src/MdtRelativeTubeT0.h
new file mode 100644
index 0000000000000000000000000000000000000000..9e5c713cf47d454c399e92f1886b54e06a59a206
--- /dev/null
+++ b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibT0/src/MdtRelativeTubeT0.h
@@ -0,0 +1,40 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#ifndef MuonCalib_MdtRelativeTubeT0_h
+#define MuonCalib_MdtRelativeTubeT0_h
+
+#include <map>
+
+namespace MuonCalib {
+
+class MdtCalibHitBase;
+class MuonFixedId;
+
+
+class MdtRelativeTubeT0
+ {
+ public:
+ enum TubeGroup {CHAMBER, MULTILAYER, LAYER, MEZZ_CARD, UNKNOWN};
+
+ inline MdtRelativeTubeT0() {}
+ inline ~MdtRelativeTubeT0() {}
+
+ void AddHit(const MdtCalibHitBase * hit);
+
+ double GetRelativeOffset(const MuonFixedId & id, TubeGroup grp) const;
+
+ private:
+
+ std::map<MuonFixedId, double> m_tube_t0;
+
+ mutable std::map<TubeGroup, std::map<MuonFixedId, double> > m_relative_offset;
+
+ inline void calculate_relative_t0s(const TubeGroup &grp) const;
+ };
+
+
+}//namespace MuonCalib
+
+#endif
diff --git a/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibT0/src/T0CalibrationClassic.cxx b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibT0/src/T0CalibrationClassic.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..c1ff600b9456bad09db2ff0fcad4008fb8b4fd26
--- /dev/null
+++ b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibT0/src/T0CalibrationClassic.cxx
@@ -0,0 +1,834 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#include "MdtCalibT0/T0CalibrationClassic.h"
+#include "MdtCalibT0/T0CalibrationOutput.h"
+#include "MuonCalibEventBase/MuonCalibSegment.h"
+#include "MuonCalibIdentifier/MuonFixedId.h"
+
+#include "MuonCalibStl/ToString.h"
+#include "MuonCalibStl/DeleteObject.h"
+
+#include <algorithm>
+#include <iostream>
+#include <vector>
+#include <string>
+
+#include "TH1.h"
+#include "TF1.h"
+#include "TH1F.h"
+#include "TH2F.h"
+#include "TFile.h"
+#include "TMinuit.h"
+
+namespace MuonCalib {
+
+//function to be fitted to time spectrum
+
+inline Double_t TimeSpectrum_func(Double_t *xx , Double_t *par)
+ {
+ Double_t &x(xx[0]);
+ return par[0]+(par[1]*(1+par[2]*exp(-(x-par[4])/par[3])))/((1+exp((-x+par[4])/par[6]))*(1+exp((x-par[5])/par[7])));
+ }
+
+ T0CalibrationClassic::T0CalibrationClassic( std::string name, const T0ClassicSettings* settings ) :
+ IMdtCalibration(name), m_settings(settings), m_converged(false), m_name(name),m_result(0), m_delete_settings(false)
+ {
+ if (!m_settings){
+ double * params=new double[8]; // warning: this would give a memory leak
+ params[0] = 0. ; // noise level
+ params[1] = 6.5 ;
+ params[2] = 6.5 ;
+ params[3] = 155. ;
+ params[4] = 280. ; //t0
+ params[5] = 980. ; //tmax
+ params[6] = 9. ; //t0 slope
+ params[7] = 11.5 ; //tmax slope
+ m_settings=new T0ClassicSettings(0.,300.,100,-100.,900.,1000,1,1000,1,8,params,10.,4);
+ m_delete_settings=true;
+ }
+ if(m_settings->printLevel() <= 3)
+ std::cout<<"T0CalibrationClassic::T0CalibrationClassic"<<m_name<<" "<<name<<std::endl;
+ if(m_settings->printLevel() <= 2) m_settings->print();
+ m_currentItnum=1;
+ std::string HistoFileName="T0Classic_"+m_name+".root";
+ if(m_settings->printLevel() <= 3)
+ std::cout<<"T0CalibrationClassic::T0CalibrationClassic"<<m_name
+ <<" "<<name<<" "<<HistoFileName<<std::endl;
+ p_file = new TFile(HistoFileName.c_str(),"recreate");
+ m_regiondir = p_file->mkdir(m_name.c_str());
+ }
+
+ T0CalibrationClassic::~T0CalibrationClassic()
+ {
+ if(m_settings->printLevel() <= 3) std::cout << "T0CalibrationClassic::~T0CalibrationClassic()" << std::endl;
+ p_file->Write();
+ p_file->Close();
+ delete p_file;
+ std::for_each(p_histos.begin(),p_histos.end(),DeleteObject());
+ if(m_delete_settings)
+ delete m_settings;
+ }
+
+ bool T0CalibrationClassic::handleSegment( MuonCalibSegment& seg )
+ {
+ for(std::vector<MdtCalibHitBase*>::iterator it =seg.mdtHOTBegin() ;
+ it!=seg.mdtHOTEnd();++it)
+ {
+ //
+ // M.I. 16 October 2007
+ // ATTENTION ATTENTION Including cut on DistanceToReadOut
+ // ONLY FOR P1 data
+ //
+ float distanceToRO = (*it)->distanceToReadout() ;
+
+ bool ROside = distanceToRO<130000. ; // this means that there is no selection along the tube
+ // bool ROside = distanceToRO<1300. ;
+// bool HVside = distanceToRO>2200. ;
+ // bool ChamberChop = distanceToRO>100.&&distanceToRO<600. ;
+ // bool ChamberChop = distanceToRO>600.&&distanceToRO<1100. ;
+ // bool ChamberChop = distanceToRO>1100.&&distanceToRO<1600. ;
+ // bool ChamberChop = distanceToRO>1900.&&distanceToRO<2400. ;
+ // bool ChamberChop = distanceToRO>2400.&&distanceToRO<2900. ;
+ // bool ChamberChop = distanceToRO>2900.&&distanceToRO<3400. ;
+
+ if (ROside) {
+
+ MuonFixedId id=(*it)->identify();
+
+ // get the T0 originally subtracted for this hit
+ int nML=id.mdtMultilayer();
+ int nL=id.mdtTubeLayer();
+ int nT=id.mdtTube();
+ // std::cout<<"Accessing Single Tube Calib info found for ML="<<nML<<" L="<<nL<<" T="<<nT<<" mezzanine="<<id.mdtMezzanine()<<std::endl;
+ const MdtTubeFitContainer::SingleTubeCalib * stc=m_result->getCalib(nML-1,nL-1,nT-1);
+// double oldT0=0;
+ if(!stc)
+ {
+ std::cout<<"no Single Tube Calib info found for ML="<<nML<<" L="<<nL<<" T="<<nT<<std::endl;
+ std::cout<<"container size "<<m_result->size()<<std::endl;
+ std::cout<<"container nML "<<m_result->numMultilayers()<<std::endl;
+ std::cout<<"container nL "<<m_result->numLayers()<<std::endl;
+ std::cout<<"container nT "<<m_result->numTubes()<<std::endl;
+ }
+
+ // get histos
+ T0ClassicHistos* histos = getHistos(id.getIdInt());
+
+ // fill histos
+ // std::cout<<"T0Classic: drifttime "<<(*it)->driftTime()<<" oldt0 "<<oldT0<<" tdc counts "<<(*it)->tdcCount()<<" propagation "<<(*it)->propagationTime()<<" slewing "<<(*it)->slewingTime()<<" flight "<<(*it)->timeOfFlight()<<std::endl;
+ // float ttime=(((*it)->tdcCount())*(25./32.))-((*it)->propagationTime())-((*it)->slewingTime())-((*it)->timeOfFlight());
+
+
+ //Luca 11/6/2007 - Tolgo il tempo di volo per i cosmici
+ // float ttime=(((*it)->tdcCount())*(25./32.))-((*it)->propagationTime())-((*it)->timeOfFlight());
+
+ //M.I. 22/6/2007 - Tolgo il tempo di propagaz per i cosmici sett.5
+ //F.R. Always use the drift time, and not the tdc-count. It will be filled correctly bu the CalibNtupleAnalysysAlg
+ //float ttime=((*it)->tdcCount())*(25./32.);
+
+ float ttime=((*it)->driftTime());
+ // histos->time->Fill((*it)->tdcCount()+oldT0);
+
+ histos->time->Fill(ttime);
+ histos->adc->Fill((*it)->adcCount());
+ //histos.adc_vs_time->Fill((*it)->driftTime(),(*it)->adcCount);
+ // book an additional dummy set of histograms to fit the global T0
+ T0ClassicHistos* histosAll = getHistos(0);
+ histosAll->time->Fill(ttime);
+ histosAll->adc->Fill((*it)->adcCount());
+
+ // M.I. Jun20-07 ---- Adding MultiLayer Histos :
+ T0ClassicHistos* histosML = getHistos(nML);
+ histosML->time->Fill(ttime);
+ histosML->adc->Fill((*it)->adcCount());
+ // M.I. Jan1507 ---- Adding Mezzanine Histos :
+ T0ClassicHistos* histosMezz = getHistos(id.mdtMezzanine());
+ histosMezz->time->Fill(ttime);
+ histosMezz->adc->Fill((*it)->adcCount());
+
+ // M.I. 16Oct07 ---- Adding "SerialGas" Histos :
+ int serialGas = nML*10+(nT-1)%3+1 ;
+ T0ClassicHistos* histosSerialGas = getHistos(serialGas);
+ histosSerialGas->time->Fill(ttime);
+ histosSerialGas->adc->Fill((*it)->adcCount());
+ }
+ }
+ return true;
+ }
+
+ bool T0CalibrationClassic::analyse()
+ {
+ if(m_settings->printLevel() <= 3)
+ std::cout << "T0CalibrationClassic::analyse iteration "<<m_currentItnum
+ << std::endl;
+
+
+ for(std::vector<T0ClassicHistos*>::iterator it =p_histos.begin() ;
+ it!=p_histos.end();++it)
+ // loop over p_histos histograms
+ {
+ if(m_settings->fitTime()) {
+ MdtTubeFitContainer::SingleTubeFit full;
+ MdtTubeFitContainer::SingleTubeCalib st;
+ int idtube=(*it)->id ;
+ // if((*it)->id!=0 && (int)(idtube/100000000)!=9 && idtube!=1 && idtube!=2 ) {
+ if((*it)->id!=0 && (int)(idtube/100000000)!=9 ) {
+ // doTimeFit(*it,full,st);
+ doTimeFit(*it,full,st);
+ doAdcFit(*it,full,st);
+ MuonFixedId fId((*it)->id);
+ int nML=fId.mdtMultilayer();
+ int nL=fId.mdtTubeLayer();
+ int nT=fId.mdtTube();
+
+ bool setInfo=m_result->setCalib(nML-1,nL-1,nT-1,st);
+ if(!setInfo)
+ std::cout<<"T0CalibrationClassic::PROBLEM! could not set SingleTubeCalib info "<<std::endl;
+ setInfo=m_result->setFit(nML-1,nL-1,nT-1,full);
+ if(!setInfo)
+ std::cout<<"T0CalibrationClassic::PROBLEM! could not set SingleTubeFullInfo info "<<std::endl;
+
+ // DA TOGLIERE !!!!
+ //
+ // int headid,lowrun,uprun;
+ // int ok=m_db->getTubeHead(headid,lowrun,uprun);
+ // std::cout<<" headid,lowrun,uprun " << headid <<" "<<lowrun<<" "<<uprun <<std::endl;
+ }
+ }
+ }
+
+ m_currentItnum++;
+ p_file->Write();
+ return true;
+ }
+
+ const IMdtCalibrationOutput* T0CalibrationClassic::analyseSegments( const std::vector<MuonCalibSegment*> & segs )
+ {
+ for(unsigned int i=0; i<segs.size(); i++) handleSegment(*segs[i]);
+ analyse();
+ return getResults();
+ }
+
+ void T0CalibrationClassic::doTimeFit(T0ClassicHistos * T0h, MdtTubeFitContainer::SingleTubeFit & fi, MdtTubeFitContainer::SingleTubeCalib & stc)
+ {
+ TMinuit *gMinuit = new TMinuit();
+ // THIS SETTINGS VARIABLE HAS TO BE IMPLEMENTED (as of March 7, 2007) :
+ // int fitMezz = m_settings->fitMezzanine();
+ // int fitMezz(123) ;
+ int fitMezz(1) ;
+
+ int np=m_settings->numParams();
+ double * pfit=new double[np];
+ double * errfit=new double[np];
+ Double_t **matrix = new Double_t*[np];
+ for(int i=0; i<np; i++)
+ {
+ matrix[i]=new Double_t[np];
+ }
+ double * pdefault=m_settings->params();
+ double chi2;
+ int ndof;
+
+ MuonFixedId fId(T0h->id);
+
+ int isMultilayer(0) ;
+
+ if (T0h->id==1 || T0h->id==2 || (T0h->id>10&&T0h->id<=23)) isMultilayer=1 ;
+
+ if (( fId.isValid() && fId.is_mdt()) || isMultilayer ) {
+
+ std::cout<<" STARTING doTimeFit "<< std::endl;
+ TH1 * h(NULL) ;
+ if (!fitMezz || isMultilayer) {
+ std::cout<< " DEBUGGGG : " << T0h->id << std::endl ;
+ h=T0h->time;
+ }
+
+ if (fitMezz==123 && !isMultilayer) { // FIT Mixed
+ h=T0h->time;
+ if (h->GetEntries()<m_settings->entries()) {
+ int hIdMezz = fId.mdtMezzanine() ;
+ ToString ts;
+ std::string HistoId(std::string("time_mezz_") + ts((hIdMezz)%(900000000)));
+ TH1F * timeHis=(TH1F*) m_regiondir->Get(HistoId.c_str());
+ if(!timeHis) {
+ delete [] pfit ;
+ delete [] errfit ;
+ for(int i=0; i<np; i++)
+ delete [] matrix[i];
+ delete [] matrix;
+ return ;
+ }
+
+ T0ClassicHistos* histosMezz = getHistos(fId.mdtMezzanine());
+ h= histosMezz->time;
+ if (h->GetEntries()<m_settings->entries()) {
+ int nML=fId.mdtMultilayer();
+ T0ClassicHistos* histosML = getHistos(nML);
+ h= histosML->time;
+ }
+ }
+ }
+
+ if (fitMezz==1 && !isMultilayer) { // FIT MEZZANINE
+ int hIdMezz = fId.mdtMezzanine() ;
+
+ ToString ts;
+ std::string HistoId(std::string("time_mezz_") + ts((hIdMezz)%(900000000)));
+ if (m_settings->printLevel() <= 2) {
+ std::cout<<" doTimeFit HistogramId : "<< T0h->id << std::endl;
+ std::cout<<" doTimeFit Histogram : "<< HistoId << std::endl;
+ }
+ TH1F * timeHis=(TH1F*) m_regiondir->Get(HistoId.c_str());
+ if(!timeHis) {
+ delete [] pfit ;
+ delete [] errfit ;
+ for(int i=0; i<np; i++)
+ delete [] matrix[i];
+ delete [] matrix;
+ return ;
+ }
+
+ T0ClassicHistos* histosMezz = getHistos(fId.mdtMezzanine());
+ h= histosMezz->time;
+ }
+
+ if (fitMezz==2 && !isMultilayer) { // FIT MULTILAYER
+ int nML=fId.mdtMultilayer();
+ T0ClassicHistos* histosML = getHistos(nML);
+ h= histosML->time;
+ }
+
+
+ Stat_t entries=h->GetEntries();
+ fi.statistics=(int)entries;
+
+ if(m_settings->printLevel() <= 1) std::cout<<" histogram "<<h->GetName()
+ <<" "<<h->GetTitle()
+ <<" entries="<<h->GetEntries()
+ <<" min entries="<<m_settings->entries()
+ <<std::endl;
+
+ // CHECK whether the Selected Histogram has enough entries
+ if ((int)entries > m_settings->entries()){
+
+ // CHECK whether the histogram has been already fitted
+ TF1 * FitFunction = h->GetFunction("TimeSpectrum");
+ if(FitFunction){
+// double chiquadro = FitFunction->GetChisquare();
+ for(int i=0; i<np; i++){
+ pfit[i] = FitFunction->GetParameter(i);
+ errfit[i] = FitFunction->GetParError(i) ;
+ }
+ chi2 = FitFunction->GetChisquare() ; // total chi2
+ ndof = FitFunction->GetNDF(); // number of degrees of freedom
+ }
+ else // The Selected Histo has NOT been fitted. Fit Now
+ {
+
+ TF1 * TimeSpectrum = new TF1("TimeSpectrum",
+ "[0]+([1]*(1+[2]*exp(-(x-[4])/[3])))/((1+exp((-x+[4])/[6]))*(1+exp((x-[5])/[7]))) ",
+ m_settings->minTime(),m_settings->maxTime());
+ // if ((int)entries > m_settings->entries())
+ // searchParams(h,&pfit[0],np);
+ for(int i=0;i<np;i++) {
+ pfit[i]=*(pdefault++);
+ if(m_settings->printLevel() <= 2)
+ std::cout<<"T0CalibrationClassic::doTimeFit initial parameter "
+ <<i<<"="<<pfit[i]<<std::endl;
+ }
+ // if ( !m_settings->initParamFlag() ) {
+ searchParams(h,&pfit[0],np);
+ if(m_settings->printLevel() <= 2) {
+ std::cout<<"T0CalibrationClassic::doTimeFit parameters after searchParams "<<std::endl;
+ for(int i=0;i<np;++i){std::cout << "i,pfit(i) "<<i<<" "<<pfit[i]<<std::endl;}}
+ std::cout <<" Sto per mettere i limiti"<<std::endl;
+ TimeSpectrum->SetParameters(pfit);
+ TimeSpectrum->SetParLimits(0,0.,5.);
+ TimeSpectrum->SetParLimits(1,0.,1000.);
+ TimeSpectrum->SetParLimits(2,0.,40.);
+ TimeSpectrum->SetParLimits(3,50.,400.);
+ TimeSpectrum->SetParLimits(4,0.,1000.);
+ // 5 parameters fit
+ TimeSpectrum->SetParLimits(5,pfit[5],pfit[5]);
+ TimeSpectrum->SetParLimits(6,pfit[6],pfit[6]);
+ TimeSpectrum->SetParLimits(7,pfit[7],pfit[7]);
+ h->Fit("TimeSpectrum","QLB");
+ // 6 parameters fit
+ TimeSpectrum->SetParLimits(5,500.,2000.);
+ TimeSpectrum->SetParLimits(6,pfit[6],pfit[6]);
+ TimeSpectrum->SetParLimits(7,pfit[7],pfit[7]);
+ h->Fit("TimeSpectrum","QLB");
+ // 7 parameters fit
+ TimeSpectrum->SetParLimits(6,4.,30.);
+ h->Fit("TimeSpectrum","QLB");
+ // final 8 parameters fit
+ TimeSpectrum->SetParLimits(6,4.,30.);
+ TimeSpectrum->SetParLimits(7,4.,30.);
+ double xmin = h->GetBinLowEdge(1);
+ double xmax = pfit[5]+250.;
+ h->Fit("TimeSpectrum","QLB","",xmin,xmax);
+
+ gMinuit->mnemat(&matrix[0][0],np) ;
+ for(int i=0; i<np; i++){
+ pfit[i] = TimeSpectrum->GetParameter(i);
+ errfit[i] = TimeSpectrum->GetParError(i) ;
+ // std::cout << " cov(ii) = " << matrix[i][i] << std::endl ;
+ }
+ chi2 = TimeSpectrum->GetChisquare() ; // total chi2
+ ndof = TimeSpectrum->GetNDF(); // number of degrees of freedom
+ }
+ // THE NEW HISTOGRAM HAS BEEN FITTED
+
+
+ if(m_settings->printLevel() <= 1)
+ std::cout<<" fit results chi2/ndof="<<chi2/ndof<<" T0="<<pfit[4]<<" err="<<errfit[4]<<std::endl;
+ if(chi2/ndof < m_settings->chi2max()) {
+ stc.statusCode=0 ;// success
+ } else {
+ stc.statusCode=3 ; // bad chi2
+ }
+ stc.t0=pfit[4];
+ fi.chi2Tdc=chi2/ndof;
+ for (int i=0 ; i<np; i++) fi.par[i]=pfit[i] ;
+
+ // NOW we get rid of the covariance matrix
+ /*******************************************************
+ int jj=0;
+ for (int i=0; i<8; i++) {
+ for (int k=0; k<i+1; k++) {
+ std::cout << " i k matrix[i][k] = "<<i<<" "<<k<<" "<< matrix[i][k] <<std::endl;
+ fi.cov[jj++] = matrix[k][i] ;
+ }
+ }
+ *******************************************************/
+ // NOW we set the the first 8 values of fi.cov to errfit
+ for (int i=0 ; i<np; i++) fi.cov[i]=errfit[i] ;
+
+ //
+ // Try 4 parameters fit now
+ // NOW COMMENTED BUT SHOULD BE TRIED !!!! (as it was in calib )
+ /*****************************************************
+ double lowt;
+ double hight;
+ double * pfd=new double[4];
+ pfd[0]=pfit[0];
+ pfd[1]=pfit[1];
+ pfd[2]=pfit[4];
+ pfd[3]=pfit[7];
+ lowt=pfit[4]-2.*pfit[6];
+ hight=pfit[4]+50.;
+ TF1 * FermiDirac = new TF1("FermiDirac",
+ "[0]+([1]/(1+exp(-(x-[3])/[2])))",
+ lowt,hight);
+ // FermiDirac->SetParameters(pfd);
+ //h->Fit("FermiDirac","LBV");
+ *****************************************************/
+
+ } else {
+ stc.statusCode=2; // too few entries
+ }
+ }
+ delete [] pfit;
+ delete [] errfit;
+ for(int i=0; i<np; i++)
+ delete [] matrix[i];
+ delete [] matrix;
+
+ std::cout<<" ENDING doTimeFit "<< std::endl;
+
+ }
+
+ void T0CalibrationClassic::doAdcFit(T0ClassicHistos * T0h, MdtTubeFitContainer::SingleTubeFit & fi, MdtTubeFitContainer::SingleTubeCalib & stc)
+ {
+ // THIS SETTINGS VARIABLE HAS TO BE IMPLEMENTED (as of March 7, 2007) :
+ // int fitMezz = m_settings->fitMezzanine();
+ // int fitMezz(123) ;
+ int fitMezz(1) ;
+
+ double chi2(0) ;
+ int ndof(0) ;
+ Double_t par[4];
+ Double_t errpar[4];
+
+ for (int ii=0; ii<4; ii++) {
+ par[ii] = 0.0 ;
+ errpar[ii] = 0.0 ;
+ }
+
+ MuonFixedId fId(T0h->id);
+ if (fId.isValid() && fId.is_mdt()) {
+ std::cout<<" doAdcFit : checking Single tube entries"<< std::endl;
+ double adcThreshold = 50. ;
+ TH1 * hcheck ;
+ hcheck = T0h->adc;
+ int nhits = static_cast<int>(hcheck->GetEntries()) ;
+
+ int adcBinThreshold = (int)((adcThreshold-hcheck->GetBinLowEdge(1))/(hcheck->GetBinWidth(1))+1) ; //
+ int nhitsAboveAdcCut = (int) hcheck->Integral(adcBinThreshold,hcheck->GetNbinsX()) ;
+ std::cout<<" doAdcFit : TotHits, nhitsAboveAdcCut "<< nhits <<" "<<nhitsAboveAdcCut<< std::endl;
+
+ fi.cov[20]=nhits ;
+ fi.cov[21]=nhitsAboveAdcCut ;
+
+ std::cout<<" STARTING doAdcFit "<< std::endl;
+
+ TH1 * h ;
+ if (!fitMezz) {
+ h=T0h->adc;
+ }
+
+ if (fitMezz==123) { // FIT Mixed
+ h=T0h->adc;
+ if (h->GetEntries()<m_settings->entries()) {
+ int hIdMezz = fId.mdtMezzanine() ;
+ ToString ts;
+ std::string HistoId(std::string("time_mezz_") + ts((hIdMezz)%(900000000)));
+ // std::cout <<" DEBUGG FITTING MEZZANINE "<<fId.mdtTube()
+ // <<" "<<fId.mdtTubeLayer()<<" "<<fId.mdtMultilayer()
+ // <<" "<<fId.mdtMezzanine() <<std::endl;
+
+ TH1F * adcHis=(TH1F*) m_regiondir->Get(HistoId.c_str());
+ if(!adcHis) {
+ return ;
+ }
+
+ T0ClassicHistos* histosMezz = getHistos(fId.mdtMezzanine());
+ h= histosMezz->adc;
+ if (h->GetEntries()<m_settings->entries()) {
+ int nML=fId.mdtMultilayer();
+
+ // std::cout <<" DEBUGG FITTING MULTILAYER "<<fId.mdtTube()
+ // <<" "<<fId.mdtTubeLayer()<<" "<<nML
+ // <<" "<<fId.mdtMezzanine() <<std::endl;
+
+ T0ClassicHistos* histosML = getHistos(nML);
+ h= histosML->adc;
+ }
+ }
+ }
+
+ if (fitMezz==1) {
+ int hIdMezz = fId.mdtMezzanine() ;
+ ToString ts;
+ std::string HistoId(std::string("charge_mezz_") + ts((hIdMezz)%(900000000)));
+ if (m_settings->printLevel() <= 2) {
+ std::cout<<" doAdcFit HistogramId : "<< T0h->id << std::endl;
+ std::cout<<" doAdcFit Histogram : "<< HistoId << std::endl;
+ }
+ TH1F * adcHis=(TH1F*) m_regiondir->Get(HistoId.c_str());
+ if(!adcHis) {
+ return ;
+ }
+ T0ClassicHistos* histosMezz = getHistos(fId.mdtMezzanine());
+ h= histosMezz->adc;
+ }
+ if (fitMezz==2) { // FIT MULTILAYER
+ int nML=fId.mdtMultilayer();
+ T0ClassicHistos* histosML = getHistos(nML);
+ h= histosML->adc;
+ }
+
+ if(m_settings->printLevel() <= 1)
+ std::cout<<" histogram "<<h->GetName()<<" "<<h->GetTitle()
+ <<" entries="<<h->GetEntries()<<std::endl;
+ Stat_t entries=h->GetEntries();
+
+ // CHECK whether the Selected Histogram has enough entries
+ if ((int)entries > m_settings->entries()){
+
+ // CHECK whether the histogram has been already fitted
+ TF1 * FitFunction = h->GetFunction("AdcSpectrum");
+ if (FitFunction){
+ chi2 = FitFunction->GetChisquare();
+ ndof = FitFunction->GetNDF();
+ for (int i=0; i<4; i++) {
+ par[i] = FitFunction->GetParameter(i);
+ errpar[i] = FitFunction->GetParError(i);
+ }
+
+ }
+ else // The Selected Histo has NOT been fitted. Fit Now
+ {
+ double m=h->GetMean();
+ double r=h->GetRMS();
+ double maxval=h->GetMaximum();
+ double adcpar[4] ;
+ adcpar[0] = maxval*2. ;
+ adcpar[1] = m ;
+ adcpar[2] = r ;
+ adcpar[3] = r/3. ;
+
+ TF1 * AdcSpectrum = new TF1("AdcSpectrum"," ([0]*exp((x-[1])/[2]))/ (1.+exp((x-[1])/[3])) ",
+ m_settings->minAdc(),m_settings->maxAdc() );
+ AdcSpectrum->SetParameters(adcpar) ;
+ double fitMin = m-(3*r) ;
+ double fitMax = m+(3*r) ;
+ if (fitMin<adcThreshold) fitMin = adcThreshold ;
+ if (fitMax>300 ) fitMax = 300. ;
+ h->Fit("AdcSpectrum","Q"," ",fitMin,fitMax);
+ chi2 = AdcSpectrum->GetChisquare();
+ ndof = AdcSpectrum->GetNDF();
+ for (int i=0; i<4; i++) {
+ par[i] = AdcSpectrum->GetParameter(i);
+ errpar[i] = AdcSpectrum->GetParError(i);
+ }
+ if (m_settings->printLevel() <=1)
+ std::cout<<"chi2/ndof="<<chi2/ndof<<" "<<"Mean="<<m<<" "<<"RMS="<<r
+ <<" par 0 1 2 3 " << par[0] <<" "<<par[1]<<" "<<par[2]<< " "<<par[3]<< std::endl;
+
+ }
+ // THE NEW HISTOGRAM HAS BEEN FITTED
+
+ } else {
+ // stc.statusCode=2; // too few entries - DO NOT CHANGE statusCode coming from doTimeFit
+ }
+
+ stc.adcCal=par[1];
+ // fi.chi2Adc=chi2/ndof; // ??? esiste fi.chi2Adc ????
+ // for (int i=0 ; i<np; i++) fi.par[i]=pfit[i] ; ...dove li mettiamo ???
+ fi.adc_par[0] = par[1] ;
+ fi.adc_par[1] = par[2] ;
+ fi.adc_err[0] = errpar[1] ;
+ fi.adc_err[1] = errpar[2] ;
+ fi.adc_chi2 = chi2/ndof ;
+ }
+ }
+
+ void T0CalibrationClassic::searchParams(TH1 * h, double * p, int np)
+ {
+ int nbinsX=h->GetNbinsX();
+ double sizeX = h->GetBinWidth(1);
+ double oldSizeX=sizeX;
+ int RebinFactor = static_cast<int>(10./sizeX);
+ // extract starting values for fit params p[np] from the Time Spectrum h
+ TH1 *hnew = h->Rebin(RebinFactor,"hnew"); // creates a new histogram hnew
+ //merging 5 bins of h1 in one bin
+ std::cout << "nbinsx,sizex,rebinfactor="<<nbinsX<<" "<<sizeX<<" "<<RebinFactor<<std::endl;
+ float minDeriv(9999.);
+ int minDerivBin(0);
+ sizeX = hnew->GetBinWidth(1);
+ int newbins = hnew->GetNbinsX();
+ for(int i=0;i<np;++i){std::cout << "i,p(i) "<<i<<" "<<p[i]<<std::endl;}
+ for(int i=0; i<newbins-1; ++i) {
+ if(hnew->GetBinContent(i)-hnew->GetBinContent(i+1)
+ <minDeriv) {
+ minDeriv = hnew->GetBinContent(i)-hnew->GetBinContent(i+1) ;
+ minDerivBin = i ;
+ }
+ }
+ float t0guess = hnew->GetBinCenter(minDerivBin);
+ if (minDerivBin<newbins-1) {
+ t0guess += (hnew->GetBinCenter(minDerivBin+1)-hnew->GetBinCenter(minDerivBin))/2.;
+ }
+ std::cout << " t0guess is "<<t0guess<<std::endl;
+ //
+ // =================== Noise level search ===================================
+ //
+ float noise(0) ;
+ int numOfBins(10), numOfBinsOffset(3) ;
+ int imin, imax ;
+ if (minDerivBin>numOfBins+numOfBinsOffset) {
+ imin = minDerivBin-numOfBins-numOfBinsOffset;
+ imax = minDerivBin-numOfBinsOffset;
+ } else {
+ imin = 0 ;
+ if (minDerivBin>numOfBinsOffset) {
+ imax = minDerivBin-numOfBinsOffset;
+ } else {
+ imax = minDerivBin;
+ }
+ }
+ int icount(0) ;
+ for (int i=imin; i<=imax; ++i) {
+ noise += hnew->GetBinContent(i);
+ icount++ ;
+ }
+
+ noise = noise/(float)(icount) ;
+ std::cout << " noise is "<<noise<<std::endl;
+ //
+ // =================== Normalization =========================================
+ //
+ int t0bin = minDerivBin;
+ int ix1 = t0bin+(int)(50/sizeX);
+ int ix2 = t0bin+(int)(500/sizeX);
+ std::cout << "t0bin,ix1,ix2 "<<t0bin<<" "<<ix1<<" "<<ix2<<std::endl;
+ float P1=p[1];
+ float P2=p[2];
+ float P3=p[3];
+ std::cout <<"P1,P2,P3 start are "<<P1<<" "<<P2<<" "<<P3<<std::endl;
+ p[0]=noise;
+ p[4]=t0guess;
+ p[5]=p[4]+700;
+ p[1]=20.;
+ p[2]=10.;
+ if (0<ix1 && ix1<newbins && 0<ix2 && ix2<newbins ) {
+ float a1=hnew->GetBinCenter(ix1);
+ float a2=hnew->GetBinCenter(ix2);
+ float cont1=hnew->GetBinContent(ix1);
+ float cont2=hnew->GetBinContent(ix2);
+ if (cont1>0. && cont2>0. ){
+ float A1=exp(-(a1-t0guess)/P3);
+ float A2=exp(-(a2-t0guess)/P3);
+ // do not forget rebinning!
+ P2 = (cont1/cont2-1.)/(A1-cont1/cont2*A2);
+ P1 = cont1/(1+P2*A1);
+ P1=P1*oldSizeX/sizeX ;
+ P2=P2*oldSizeX/sizeX ;
+ std::cout << "a1,a2 "<<a1<<" "<<a2<<" cont1,cont2 "<<cont1<<" "<<cont2<<" A1,A2 "<<A1<<" "<<A2<<std::endl;
+ std::cout << " t0Guess .... P1, P2 " <<P1<<" "<<P2<<std::endl;
+ p[1]=P1;
+ p[2]=P2;
+ }
+ }
+ delete hnew;
+ return;
+ }
+
+ const IMdtCalibrationOutput* T0CalibrationClassic::getResults() const
+ {
+ return new T0CalibrationOutput(m_result);
+ }
+
+ T0ClassicHistos* T0CalibrationClassic::getHistos(unsigned int idtube)
+ {
+ ToString ts;
+ // std::cout<<"T0CalibrationClassic::getHistos "<<idtube<<std::endl;
+ std::string HistoId;
+ if((int)(idtube/100000000)==9){
+ // std::cout<<"T0CalibrationClassic::getHistos "<<idtube<<" of tipe mezz(idtube)%(900000000), mezz "<<(idtube)%(900000000)<<std::endl;
+ int mezz=(idtube)%(900000000);
+ HistoId="time_mezz_"+ts(mezz);
+ } else if(idtube==0){
+ // std::cout<<"T0CalibrationClassic::getHistos "<<idtube<<" of tipe 0"<<std::endl;
+ HistoId="time";
+ } else if(idtube==1){
+ HistoId="time_ML1";
+ } else if(idtube==2){
+ HistoId="time_ML2";
+ } else if(idtube==11){
+ HistoId="time_ML1_series1";
+ } else if(idtube==12){
+ HistoId="time_ML1_series2";
+ } else if(idtube==13){
+ HistoId="time_ML1_series3";
+ } else if(idtube==21){
+ HistoId="time_ML2_series1";
+ } else if(idtube==22){
+ HistoId="time_ML2_series2";
+ } else if(idtube==23){
+ HistoId="time_ML2_series3";
+ } else {
+ // std::string HistoId(std::string("time_") + ts(idtube));
+ MuonFixedId fId(idtube);
+ int nML=fId.mdtMultilayer();
+ int nL=fId.mdtTubeLayer();
+ int nT=fId.mdtTube();
+ int tubeid=(nML*1000)+(nL*100)+nT;
+ int nstat=fId.stationName();
+ std::string stationName=fId.stationNumberToFixedStationString(nstat);
+ int eta = fId.eta();
+ int phi = fId.phi();
+ HistoId="time_"+stationName+"_"+ts(eta)+"_"+ts(phi)+"_"+ts(tubeid);
+ // std::cout<<"T0CalibrationClassic::getHistos "<<idtube<<" of tipe tube "<<HistoId<<std::endl;
+ }
+ TH1F * timeHis=(TH1F*) m_regiondir->Get(HistoId.c_str());
+ if(!timeHis) {
+ // std::cout<<"T0CalibrationClassic::getHistos "<<idtube<<" pointer not found, booking"<<std::endl;
+ return bookHistos(idtube);
+ } else {
+ for(std::vector<T0ClassicHistos*>::iterator it =p_histos.begin() ;
+ it!=p_histos.end();++it)
+ // loop over p_histos histograms to look for the set of histos of tube idtube
+ {
+ if((*it)->time==timeHis) return (*it);
+ }
+ }
+ return new T0ClassicHistos();
+ }
+
+ T0ClassicHistos* T0CalibrationClassic::bookHistos(unsigned int idtube)
+ {
+ T0ClassicHistos * histos=new T0ClassicHistos();
+ ToString ts;
+ std::string histonametdc;
+ std::string histonameadc;
+
+ histos->id=idtube;
+ m_regiondir->cd();
+ if((int)(idtube/100000000)==9){
+ int mezz=(idtube)%(900000000);
+ histonametdc="time_mezz_"+ts(mezz);
+ histonameadc="charge_mezz_"+ts(mezz);
+ } else if(idtube==0){
+ histonametdc="time";
+ histonameadc="charge";
+ } else if(idtube==1){
+ histonametdc="time_ML1";
+ histonameadc="charge_ML1";
+ } else if(idtube==2){
+ histonametdc="time_ML2";
+ histonameadc="charge_ML2";
+ } else if(idtube==11){
+ histonametdc="time_ML1_series1";
+ histonameadc="charge_ML1_series1";
+ } else if(idtube==12){
+ histonametdc="time_ML1_series2";
+ histonameadc="charge_ML1_series2";
+ } else if(idtube==13){
+ histonametdc="time_ML1_series3";
+ histonameadc="charge_ML1_series3";
+ } else if(idtube==21){
+ histonametdc="time_ML2_series1";
+ histonameadc="charge_ML2_series1";
+ } else if(idtube==22){
+ histonametdc="time_ML2_series2";
+ histonameadc="charge_ML2_series2";
+ } else if(idtube==23){
+ histonametdc="time_ML2_series3";
+ histonameadc="charge_ML2_series3";
+ } else {
+ MuonFixedId fId(idtube);
+ int nML=fId.mdtMultilayer();
+ int nL=fId.mdtTubeLayer();
+ int nT=fId.mdtTube();
+ int nstat=fId.stationName();
+ int tubeid=(nML*1000)+(nL*100)+nT;
+ std::string stationName=fId.stationNumberToFixedStationString(nstat);
+ int eta = fId.eta();
+ int phi = fId.phi();
+ histonametdc="time_"+stationName+"_"+ts(eta)+"_"+ts(phi)+"_"+ts(tubeid);
+ histonameadc="charge_"+stationName+"_"+ts(eta)+"_"+ts(phi)+"_"+ts(tubeid);
+ }
+
+ histos->time=new TH1F(histonametdc.c_str(),"Drift Time",m_settings->binTime(),m_settings->minTime(),m_settings->maxTime());
+
+ histos->adc= new TH1F(histonameadc.c_str(),"ADC",m_settings->binAdc(),m_settings->minAdc(),m_settings->maxAdc());
+
+ p_histos.push_back(histos);
+ return histos;
+ }
+
+ void T0CalibrationClassic::setInput( const IMdtCalibrationOutput* calib_in )
+ {
+
+ // This method is called both by the event loop and by the tool.
+ // Only the call from the tool is relevant for this implementation
+ // and should be performed only once.
+
+ if(m_result) return;
+
+ const T0CalibrationOutput* t0Input = dynamic_cast<const T0CalibrationOutput*>(calib_in);
+ m_result = t0Input->t0s();
+ m_result->setImplementation("T0CalibrationClassic");
+ }
+
+ bool T0CalibrationClassic::converged() const
+ {
+ return m_converged;
+ }
+
+}
diff --git a/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibT0/src/T0CalibrationMT.cxx b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibT0/src/T0CalibrationMT.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..f88dab04e8edac394d6fa20056059b8d3a8ab764
--- /dev/null
+++ b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibT0/src/T0CalibrationMT.cxx
@@ -0,0 +1,343 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#include "MdtCalibT0/T0CalibrationMT.h"
+#include "MdtCalibT0/T0CalibrationOutput.h"
+#include "MdtRelativeTubeT0.h"
+#include "MuonCalibEventBase/MuonCalibSegment.h"
+#include "MuonCalibIdentifier/MuonFixedId.h"
+#include "MdtCalibT0/HistogramId.h"
+
+#include "MdtCalibData/MdtTubeFitContainer.h"
+#include "MdtCalibT0/T0MTSettings.h"
+#include "MdtCalibT0/T0MTHistos.h"
+#include "MdtCalibT0/ADCMTHistos.h"
+
+#include "MuonCalibStl/ToString.h"
+#include "MuonCalibStl/DeleteObject.h"
+
+#include <algorithm>
+#include <iostream>
+#include <vector>
+#include <string>
+
+#include "TH1.h"
+#include "TF1.h"
+#include "TH1F.h"
+#include "TH2F.h"
+#include "TFile.h"
+
+namespace MuonCalib {
+
+
+
+T0CalibrationMT::T0CalibrationMT( std::string name, const T0MTSettings* settings , const std::vector<int> &sort_by, const std::vector<int> &adc_sort_by) :
+ IMdtCalibration(name), m_settings(settings), m_converged(false), m_name(name),m_currentItnum(0), m_sort_by(sort_by), m_adc_sort_by(adc_sort_by), m_delete_settings(false)
+{
+ if (!m_settings){
+ m_settings=new T0MTSettings();
+ m_delete_settings = true;
+ }
+
+ std::string HistoFileName="T0MT_"+m_name+".root";
+ p_file = new TFile(HistoFileName.c_str(),"recreate");
+ m_regiondir = p_file->mkdir(m_name.c_str());
+
+ p_histos.resize(sort_by.size());
+ p_adc_histos.resize(adc_sort_by.size());
+
+ m_tube_ids.resize(sort_by.size());
+ m_adc_tube_ids.resize(adc_sort_by.size());
+}
+
+T0CalibrationMT::~T0CalibrationMT()
+{
+ p_file->Write();
+ p_file->Close();
+ delete p_file;
+ for(unsigned int i=0; i<p_histos.size(); i++)
+ for(std::map<HistogramId, T0MTHistos*> :: iterator it = p_histos[i].begin(); it!= p_histos[i].end(); it++)
+ delete it->second;
+if (m_delete_settings)
+ delete m_settings;
+}
+
+bool T0CalibrationMT::handleSegment( MuonCalibSegment& seg )
+{
+ for(std::vector<MdtCalibHitBase*>::iterator it =seg.mdtHOTBegin() ;
+ it!=seg.mdtHOTEnd();++it)
+ {
+ MuonFixedId id=(*it)->identify();
+ m_nhits_per_tube[id.getIdInt()]++;
+ // get the T0 originally subtracted for this hit
+ int nML=id.mdtMultilayer();
+ int nL=id.mdtTubeLayer();
+ int nT=id.mdtTube();
+ const MdtTubeFitContainer::SingleTubeCalib * stc(NULL);
+ NtupleStationId sid(id);
+ sid.SetMultilayer(0);
+ std::map<NtupleStationId, MdtTubeFitContainer *> :: const_iterator res_it(m_result.find(sid));
+ if(res_it!=m_result.end())
+ stc=res_it->second->getCalib(nML-1,nL-1,nT-1);
+ double oldT0=0;
+ if(stc)
+ oldT0 = stc->t0;
+ else
+ {
+ std::cout<<"no Single Tube Calib info found for ML="<<nML<<" L="<<nL<<" T="<<nT<<std::endl;
+// std::cout<<"container size "<<m_result->size()<<std::endl;
+// std::cout<<"container nML "<<m_result->numMultilayers()<<std::endl;
+// std::cout<<"container nL "<<m_result->numLayers()<<std::endl;
+// std::cout<<"container nT "<<m_result->numTubes()<<std::endl;
+ }
+
+ // get histos
+ for(unsigned int i=0; i<m_sort_by.size(); i++)
+ {
+ T0MTHistos* histos = getHistos(id, i);
+ histos->GetTSpec()->Fill((*it)->driftTime() + oldT0 + (*it)->tubeT0());
+ }
+ for(unsigned int i=0; i<m_adc_sort_by.size(); i++)
+ {
+ ADCMTHistos * adc_histos = getADCHistos(id, i);
+ adc_histos->GetADCSpec()->Fill((*it)->adcCount());
+ }
+ //relative t0 offsets
+ m_rel_tube_t0s[sid].AddHit(*it);
+ }
+ return true;
+}
+
+bool T0CalibrationMT::analyse()
+ {
+ std::map<int, MdtTubeFitContainer::SingleTubeFit> full;
+ std::map<int, MdtTubeFitContainer::SingleTubeCalib> st;
+ std::map<int, std::string> fit_by;
+ if(m_settings->FitTime())
+ for(unsigned int i=0; i<m_sort_by.size(); i++)
+ {
+ analyse_tdc(i, full, st, fit_by);
+ }
+ for(unsigned int i=0; i<m_adc_sort_by.size(); i++)
+ {
+ analyse_adc(i, full, st);
+ }
+
+
+ for(std::map<int, MdtTubeFitContainer::SingleTubeFit>::iterator it=full.begin(); it!=full.end(); it++)
+ {
+ if(it->first == 0) continue;
+ MuonFixedId fId(it->first);
+ NtupleStationId sid(fId);
+ sid.SetMultilayer(0);
+ MdtTubeFitContainer::SingleTubeFit &fi(it->second);
+ MdtTubeFitContainer::SingleTubeCalib &stc(st[it->first]);
+ fi.group_by=fit_by[it->first];
+ int nML=fId.mdtMultilayer();
+ int nL=fId.mdtTubeLayer();
+ int nT=fId.mdtTube();
+ bool setInfo=m_result[sid]->setCalib(nML-1,nL-1,nT-1,stc);
+ if(!setInfo)
+ std::cout<<"T0CalibrationMT::PROBLEM! could not set SingleTubeCalib info "<<std::endl;
+ fi.n_hits=m_nhits_per_tube[fId.getIdInt()];
+ fi.n_hits_above_adc_cut=m_nhits_per_tube[fId.getIdInt()];
+ setInfo=m_result[sid]->setFit(nML-1,nL-1,nT-1,fi);
+ if(!setInfo)
+ std::cout<<"T0CalibrationMT::PROBLEM! could not set SingleTubeFit info "<<std::endl;
+ }
+
+ m_currentItnum++;
+ p_file->Write();
+ return true;
+ }
+
+bool T0CalibrationMT::analyse_tdc(const int & nr, std::map<int, MdtTubeFitContainer::SingleTubeFit> & full, std::map<int, MdtTubeFitContainer::SingleTubeCalib> & st, std::map<int, std::string> & fit_by_map)
+ {
+ if(m_settings->VerboseLevel() > 1)
+ std::cout << "T0CalibrationMT::analyse iteration "<<m_currentItnum << std::endl;
+ std::string fit_by("UNKNOWN");
+ switch(m_sort_by[nr])
+ {
+ case HistogramId::TUBE:
+ fit_by="TUBE";
+ break;
+ case HistogramId::LAYER:
+ fit_by="LAYER";
+ break;
+ case HistogramId::MULTILAYER:
+ fit_by="MULTILAYER";
+ break;
+ case HistogramId::CHAMBER:
+ fit_by="CHAMBER";
+ break;
+ case HistogramId::MEZZ_CARD:
+ fit_by="MEZZ_CARD";
+ break;
+ }
+
+ for(std::map<HistogramId, T0MTHistos*>::iterator it =p_histos[nr].begin() ;it!=p_histos[nr].end();++it)
+// loop over p_histos histograms
+ {
+ doTimeFit(it->second, m_tube_ids[nr][it->first], full, st, fit_by_map, fit_by);
+
+ }
+ return true;
+ }
+
+
+bool T0CalibrationMT::analyse_adc(const int & nr, std::map<int, MdtTubeFitContainer::SingleTubeFit> & full, std::map<int, MdtTubeFitContainer::SingleTubeCalib> & st)
+ {
+ for(std::map<HistogramId, ADCMTHistos*>::iterator it =p_adc_histos[nr].begin() ;it!=p_adc_histos[nr].end();++it)
+ if(m_settings->FitADC())
+ {
+ doAdcFit(it->second, m_adc_tube_ids[nr][it->first], full, st);
+ }
+ return true;
+ }
+
+const IMdtCalibrationOutput* T0CalibrationMT::analyseSegments( const std::vector<MuonCalibSegment*> & segs )
+ {
+ for(unsigned int i=0; i<segs.size(); i++)
+ handleSegment(*segs[i]);
+ analyse();
+ return getResults();
+ }
+
+void T0CalibrationMT::doTimeFit(T0MTHistos * T0h, const std::set<MuonFixedId> & tube_ids, std::map<int, MdtTubeFitContainer::SingleTubeFit> & fim, std::map<int, MdtTubeFitContainer::SingleTubeCalib> & stcm, std::map<int, std::string> & fit_by_map, const std::string & fit_by)
+{
+ TDirectory *cwd=gDirectory;
+ m_regiondir->cd();
+//do t0 fit
+ MdtRelativeTubeT0::TubeGroup theGroup(MdtRelativeTubeT0::UNKNOWN);
+ if(fit_by=="CHAMBER") theGroup=MdtRelativeTubeT0::CHAMBER;
+ else if(fit_by=="MULTILAYER") theGroup=MdtRelativeTubeT0::MULTILAYER;
+ else if(fit_by=="LAYER") theGroup=MdtRelativeTubeT0::LAYER;
+ else if(fit_by=="MEZZ_CARD") theGroup=MdtRelativeTubeT0::MEZZ_CARD;
+ if(T0h->FitT0() && m_settings->MinEntriesTime()<=T0h->GetTSpec()->GetEntries())
+ {
+ const TF1 *fun(T0h->GetT0Function());
+ for(std::set<MuonFixedId>::const_iterator it=tube_ids.begin(); it!= tube_ids.end(); it++)
+ {
+ if(it->getIdInt()==0) continue;
+ NtupleStationId stId(*it);
+ stId.SetMultilayer(0);
+ double rel_t0(0.0);
+ if(m_settings->T0Settings()->CorrectRelT0s())
+ rel_t0 = m_rel_tube_t0s[stId].GetRelativeOffset(*it, theGroup);
+ std::cout<<"HHhhHH "<<fit_by<<" "<<rel_t0<<std::endl;
+ MdtTubeFitContainer::SingleTubeFit &fi(fim[it->getIdInt()]);
+ MdtTubeFitContainer::SingleTubeCalib &stc(stcm[it->getIdInt()]);
+ //store parameters
+ fi.statistics = static_cast<int>(T0h->GetTSpec()->GetEntries());
+ fi.chi2Tdc=T0h->T0Chi2();
+ fi.par[0]=fun->GetParameter(T0MTHistos :: T0_PAR_NR_BACK);
+ fi.cov[0]=fun->GetParError(T0MTHistos :: T0_PAR_NR_BACK);
+ fi.par[4]=fun->GetParameter(T0MTHistos :: T0_PAR_NR_T0) + rel_t0;
+ fi.cov[4]=fun->GetParError(T0MTHistos :: T0_PAR_NR_T0);
+ fi.par[6]=fun->GetParameter(T0MTHistos :: T0_PAR_NR_T);
+ fi.cov[6]=fun->GetParError(T0MTHistos :: T0_PAR_NR_T);
+ stc.t0=fun->GetParameter(T0MTHistos :: T0_PAR_NR_T0) + rel_t0;
+ stc.statusCode=T0h->StatusCode();
+ m_converged=true;
+ fit_by_map[it->getIdInt()]=fit_by;
+ }
+ }
+ if(T0h->FitTmax()&& m_settings->MinEntriesTime()<=T0h->GetTSpec()->GetEntries())
+ {
+ const TF1 *fun(T0h->GetTMaxFunction());
+ //store parameters
+ for(std::set<MuonFixedId>::const_iterator it=tube_ids.begin(); it!= tube_ids.end(); it++)
+ {
+ if(it->getIdInt()==0) continue;
+ MdtTubeFitContainer::SingleTubeFit &fi(fim[it->getIdInt()]);
+// MdtTubeFitContainer::SingleTubeCalib &stc(stcm[*it]);
+ fi.par[5]=fun->GetParameter(T0MTHistos :: TMAX_PAR_NR_TMAX);
+ fi.cov[5]=fun->GetParError(T0MTHistos :: TMAX_PAR_NR_TMAX);
+ fi.chi2TdcEnd=fun->GetChisquare()/fun->GetNDF();
+ }
+ }
+ cwd->cd();
+ }
+
+void T0CalibrationMT::doAdcFit(ADCMTHistos * T0h, const std::set<MuonFixedId> & tube_ids, std::map<int, MdtTubeFitContainer::SingleTubeFit> & fim, std::map<int, MdtTubeFitContainer::SingleTubeCalib> & stcm)
+ {
+ if(T0h->FitAdc() && m_settings->MinEntriesADC()<=T0h->GetADCSpec()->GetEntries())
+ {
+ const TF1 *fun(T0h->GetAdcFunction());
+ if(fun==NULL) return;
+ for(std::set<MuonFixedId>::const_iterator it=tube_ids.begin(); it!= tube_ids.end(); it++)
+ {
+ if(it->getIdInt()==0) continue;
+ MdtTubeFitContainer::SingleTubeFit &fi(fim[it->getIdInt()]);
+ MdtTubeFitContainer::SingleTubeCalib &stc(stcm[it->getIdInt()]);
+
+ stc.adcCal=fun->GetParameter(1);
+ for(int i=0; (i<fun->GetNpar() && i<4); i++)
+ {
+ fi.adc_par[i]=fun->GetParameter(i);
+ fi.adc_err[i]=fun->GetParError(i);
+ }
+ fi.adc_chi2=fun->GetChisquare()/fun->GetNDF();
+ }
+ }
+ }
+
+
+const IMdtCalibrationOutput* T0CalibrationMT::getResults() const
+{
+ return new T0CalibrationOutput(m_result);
+}
+
+
+ADCMTHistos* T0CalibrationMT::getADCHistos(const MuonFixedId & idtube, unsigned int nr)
+{
+ HistogramId id;
+ id.Initialize(idtube, m_adc_sort_by[nr]);
+ if(p_adc_histos[nr][id] == NULL)
+ {
+ TDirectory *cwd=gDirectory;
+ m_regiondir->cd();
+ p_adc_histos[nr][id] = new ADCMTHistos(id.getIdInt(), m_settings, id.HistogramName().c_str());
+ cwd->cd();
+ }
+ m_adc_tube_ids[nr][id].insert(idtube);
+ return p_adc_histos[nr][id];
+}
+
+T0MTHistos* T0CalibrationMT::getHistos(const MuonFixedId & idtube, unsigned int nr) {
+ HistogramId id;
+ id.Initialize(idtube, m_sort_by[nr]);
+ if(p_histos[nr][id] == NULL)
+ {
+ TDirectory *cwd=gDirectory;
+ m_regiondir->cd();
+ p_histos[nr][id] = new T0MTHistos(id.getIdInt(), m_settings, id.HistogramName().c_str());
+ cwd->cd();
+ }
+ m_tube_ids[nr][id].insert(idtube);
+ return p_histos[nr][id];
+}
+
+void T0CalibrationMT::setInput( const IMdtCalibrationOutput* calib_in )
+{
+
+ // This method is called both by the event loop and by the tool.
+ // Only the call from the tool is relevant for this implementation
+ // and should be performed only once.
+
+ if(m_result.size()) return;
+
+ const T0CalibrationOutput* t0Input = dynamic_cast<const T0CalibrationOutput*>(calib_in);
+ if(!calib_in || !t0Input) return;
+ m_result = t0Input->GetMap();
+ for(std::map<NtupleStationId, MdtTubeFitContainer *>::iterator it=m_result.begin(); it!=m_result.end(); it++)
+ it->second->setImplementation("T0CalibrationMT");
+}
+
+bool T0CalibrationMT::converged() const
+{
+ return m_converged;
+}
+
+}
diff --git a/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibT0/src/T0MTHistos.cxx b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibT0/src/T0MTHistos.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..2d504e2836e44ed0e66f10ec399ae8ccc28b5325
--- /dev/null
+++ b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibT0/src/T0MTHistos.cxx
@@ -0,0 +1,496 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+using namespace std;
+
+#include "MdtCalibT0/T0MTHistos.h"
+#include "MdtCalibT0/MTT0PatternRecognition.h"
+#include "MdtCalibT0/MTTmaxPatternRecognition.h"
+
+//root
+#include "TLine.h"
+//root
+#include "TH1.h"
+#include "TF1.h"
+#include "TMath.h"
+#include "TDirectory.h"
+#include "TRandom.h"
+#include "list"
+#include <cmath>
+
+namespace MuonCalib {
+
+/** The fermi function to be fitted at the rising edge of the spectrum
+
+*/
+ inline Double_t mt_t0_fermi(Double_t *x , Double_t *par)
+ {
+ //more convenient parameters
+ const Double_t &t(x[0]);
+ const Double_t &t_0(par[T0MTHistos :: T0_PAR_NR_T0]),
+ &T(par[T0MTHistos :: T0_PAR_NR_T]),
+ &back(par[T0MTHistos :: T0_PAR_NR_BACK]),
+ &A(par[T0MTHistos :: T0_PAR_NR_A]);
+ //the formula
+ return (back + A/(1+exp(-(t-t_0)/T)));
+ }
+
+
+/** The fermi function to be fitted at the trailing slope of the spectrum
+
+*/
+ inline Double_t mt_tmax_fermi(Double_t *x, Double_t *par)
+ {
+ //more convenient parameters
+ Double_t &t(x[0]);
+ Double_t &t_max(par[T0MTHistos :: TMAX_PAR_NR_TMAX]),
+ &T(par[T0MTHistos :: TMAX_PAR_NR_T]),
+ &back(par[T0MTHistos :: TMAX_PAR_NR_BACK]),
+ &a(par[T0MTHistos :: TMAX_PAR_NR_A]),
+ &b(par[T0MTHistos :: TMAX_PAR_NR_B]),
+ &t_0(par[T0MTHistos :: TMAX_PAR_NR_T0]);
+ //the formula
+ return (back + (exp(a+b*(t-t_0)))/(1+exp((t-t_max)/T)));
+ }
+
+
+
+//////////////////
+// Fill.. //
+//////////////////
+
+void T0MTHistos :: FillT(double t)
+ {m_time->Fill(static_cast<Axis_t>(t));}
+
+
+//////////////////////////////////////////////////////////
+// Initialize(int id, const T0MTSettings & settings) //
+//////////////////////////////////////////////////////////
+
+void T0MTHistos :: Initialize(int id, const T0MTSettings * settings, const char * hname)
+ {
+ m_settings=settings;
+// if(m_settings->VerboseLevel()>1)
+// {
+ cout<<"T0MTHistos :: Initialize: called"<<endl;
+// }
+ char buf[100];
+ if(hname==NULL)
+ snprintf(buf, 100, "t_spec_%d", id);
+ else
+ snprintf(buf, 100, "t_spec_%s", hname);
+ std::cout<<gDirectory->GetName()<<std::endl;
+ m_time=new TH1F(buf, "", settings->NBinsTime(), settings->TimeMin(), settings->TimeMax());
+ m_id=id;
+ if(settings->DrawDebugGraphs())
+ {
+ if(m_settings->VerboseLevel()>1)
+ {
+ cout<<"T0MTHistos :: Initialize: debug directory created"<<endl;
+ }
+ TDirectory *cwd=gDirectory;
+ snprintf(buf, 100, "t0_tmax_dir_%d", id);
+ dir=gDirectory->mkdir(buf, buf);
+ cwd->cd();
+ }
+ else dir=NULL;
+ m_t0_ok=false;
+ m_tmax_ok=false;
+ }
+
+//////////////////////////////////
+// SetTSpec(int id, TH1F *spec) //
+//////////////////////////////////
+
+void T0MTHistos :: SetTSpec(int id, TH1F *spec, const T0MTSettings * settings, bool copy_spec)
+ {
+ m_settings=settings;
+ TDirectory *cwd=gDirectory;
+ if(copy_spec)
+ {
+ m_time=new TH1F(*spec);
+ }
+ else
+ {
+ m_time=spec;
+ }
+ m_id=id;
+ if(settings->DrawDebugGraphs())
+ {
+ char buffer[100];
+ snprintf(buffer, 100, "t0_tmax_dir_%d", id);
+ dir=gDirectory->mkdir(buffer, buffer);
+ cwd->cd();
+ }
+ }
+
+
+
+bool T0MTHistos :: FitT0()
+ {
+ if(m_time->GetEntries()<1000)
+ {
+ m_status_code=2;
+ return false;
+ }
+ if(!NormalFit()) return false;
+ if(m_time->GetEntries()<10000)
+ {
+ m_status_code=0;
+ return true;
+ }
+ if(m_settings->T0Settings()->ScrambleThreshold()>0 && m_chi2>m_settings->T0Settings()->ScrambleThreshold())
+ {
+ if(!T0Scramble())
+ {
+ m_status_code=3;
+ return false;
+ }
+ }
+ if(m_settings->T0Settings()->SlicingThreshold() >0 && m_chi2>m_settings->T0Settings()->SlicingThreshold())
+ {
+ TopSlicing();
+ }
+ std::cout<<m_time->GetName()<<" "<<m_chi2<<std::endl;
+ m_status_code=0;
+ return true;
+ }
+
+
+//////////////////
+// FitTmax() //
+//////////////////
+
+bool T0MTHistos :: FitTmax()
+ {
+ TDirectory *cwd=gDirectory;
+ if(dir!=NULL) dir->cd();
+ if(m_time==NULL)
+ {
+ cerr<<"T0MTHistos :: FitTmax: Class is not initialized!"<<endl;
+ m_tmax_ok=false;
+ cwd->cd();
+ return false;
+ }
+//check if t0-fit was successfull t0 is needed for tmax-pattern recognition
+ if(!m_t0_ok || m_t0_fermi==NULL)
+ {
+ cerr<<"T0MTHistos :: FitTmax for tube "<<m_id<<": No valid t0-value!"<<endl;
+ m_tmax_ok=false;
+ cwd->cd();
+ return false;
+ }
+//Create pattern Recognition Class
+ MTTmaxPatternRecognition rec;
+//perform pattern recognition
+ if(!rec.Initialize(m_time, m_t0_fermi->GetParameter(T0_PAR_NR_T0), m_settings))
+ {
+ cerr<<"T0MTHistos :: FitTmax for tube "<<m_id<<": Pattern recognition failed!"<<endl;
+ m_tmax_ok=false;
+ cwd->cd();
+ return false;
+ }
+//create function object
+// sprintf(buffer,"mt_tmax_fermi_%d", m_id);
+ char buffer[100];
+ snprintf(buffer,100, "mt_tmax_fermi");
+ if(!m_tmax_fermi){
+ m_tmax_fermi = new TF1(buffer, mt_tmax_fermi, rec.GetFitRangeMin(), rec.GetFitRangeMax(), N_TMAX_FIT_PAR);
+//set parameter names
+ m_tmax_fermi->SetParName(TMAX_PAR_NR_TMAX, "t_{max}");
+ m_tmax_fermi->SetParName(TMAX_PAR_NR_T, "T");
+ m_tmax_fermi->SetParName(TMAX_PAR_NR_BACK, "r_{b}");
+ m_tmax_fermi->SetParName(TMAX_PAR_NR_A, "a");
+ m_tmax_fermi->SetParName(TMAX_PAR_NR_B, "b");
+//set fixed values
+// cout<<"fixing r_b="<<rec.GetBackground()<<endl;
+ m_tmax_fermi->FixParameter(TMAX_PAR_NR_BACK, rec.GetBackground());
+// cout<<"fixing a="<<rec.GetA()<<endl;
+ m_tmax_fermi->FixParameter(TMAX_PAR_NR_A, rec.GetA());
+// cout<<"fixing b="<<rec.GetB()<<endl;
+ m_tmax_fermi->FixParameter(TMAX_PAR_NR_B, rec.GetB());
+ m_tmax_fermi->FixParameter(TMAX_PAR_NR_T0, m_t0_fermi->GetParameter(T0_PAR_NR_T0));
+//set start values
+// cout<<"start value for t_max="<<rec.GetEstimatedTMax()<<endl;
+ m_tmax_fermi->SetParameter(TMAX_PAR_NR_TMAX, rec.GetEstimatedTMax());
+ m_tmax_fermi->SetParameter(TMAX_PAR_NR_T, 3.0);
+ }
+//perform fit
+ if(dir!=NULL)
+ {
+ m_tmax_fermi->SetLineColor(3);
+ m_tmax_fermi->Write();
+ }
+ m_tmax_fermi->SetLineColor(4);
+ string fitopt("LR");
+ if(m_settings->VerboseLevel()==0)
+ fitopt+="Q";
+ if(m_settings->AddFitfun())
+ {
+ fitopt+="+";
+ }
+ else
+ {
+ fitopt+="0";
+ }
+ m_time->Fit(m_tmax_fermi, fitopt.c_str());
+ m_tmax_ok=true;
+ cwd->cd();
+ return true;
+ }
+
+
+
+
+//////////////////
+// FitT0() //
+//////////////////
+
+bool T0MTHistos :: NormalFit()
+ {
+ if(m_settings->VerboseLevel()>1)
+ {
+ cout<<"T0MTHistos :: FitT0(): called"<<endl;
+ }
+ TDirectory *cwd=gDirectory;
+ if(dir!=NULL) dir->cd();
+//check if class is initialized
+ if(m_time==NULL)
+ {
+ cerr<<"T0MTHistos :: FitT0: Class is not initialized!"<<endl;
+ m_t0_ok=false;
+ cwd->cd();
+ m_status_code=3;
+ return false;
+ }
+//create pattern recognition class
+ MTT0PatternRecognition rec;
+//perform pattern recognition
+ if(!rec.Initialize(m_time, m_settings))
+ {
+ m_t0_ok=false;
+ cerr<<"T0MTHistos :: FitT0 for tube "<<m_id<<": Pattern recognition failed!"<<endl;
+ cwd->cd();
+ m_status_code=3;
+ return false;
+ }
+ m_t0_ok=true;
+//create function class
+ char buffer[100];
+// sprintf(buffer,"mt_t0_fermi_%d", m_id);
+ snprintf(buffer,100, "mt_t0_fermi");
+ m_t0_fermi=new TF1(buffer, mt_t0_fermi, rec.GetFitRangeMin(), rec.GetFitRangeMax(), N_T0_FIT_PAR);
+//set parameter names
+ m_t0_fermi->SetParName(T0_PAR_NR_T0, "t_{0}");
+ m_t0_fermi->SetParName(T0_PAR_NR_T, "T");
+ m_t0_fermi->SetParName(T0_PAR_NR_BACK, "r_{b}");
+ m_t0_fermi->SetParName(T0_PAR_NR_A, "A");
+//set fixed values
+ m_t0_fermi->FixParameter(T0_PAR_NR_BACK, rec.GetBackground());
+ m_t0_fermi->SetParameter(T0_PAR_NR_A, rec.GetHeight() - rec.GetBackground());
+//set estimates as start values
+ m_t0_fermi->SetParameter(T0_PAR_NR_T0, rec.GetEstimatedT0());
+//set resonable start value for T
+ m_t0_fermi->SetParameter(T0_PAR_NR_T, 3.0);
+//perform fit - NOTE: The return value of the Fit function is not documented!
+ if(dir!=NULL)
+ {
+ TLine *ln = new TLine( rec.GetFitRangeMin(), 0, rec.GetFitRangeMin(), m_time->GetMaximum());
+ ln->Write("t0_range_min");
+ ln = new TLine( rec.GetFitRangeMax(), 0, rec.GetFitRangeMax(), m_time->GetMaximum());
+ ln->Write("t0_range_max");
+ m_t0_fermi->SetLineColor(3);
+ m_t0_fermi->Write();
+ }
+ m_t0_fermi->SetLineColor(2);
+ string fitopt("BLR");
+ if(m_settings->VerboseLevel()==0)
+ fitopt+="Q";
+ if(m_settings->AddFitfun())
+ {
+ fitopt+="+";
+ }
+ else
+ {
+ fitopt+="0";
+ }
+ m_time->Fit(m_t0_fermi, fitopt.c_str(),"", rec.GetFitRangeMin(), rec.GetFitRangeMax());
+ if(m_settings->T0Settings()->UseTopChi2())
+ TopChi2();
+ else
+ m_chi2=m_time->GetFunction("mt_t0_fermi")->GetChisquare() / m_time->GetFunction("mt_t0_fermi")->GetNDF();
+ m_t0_ok=true;
+ cwd->cd();
+ m_status_code=0;
+ return true;
+ }
+
+
+//////////////////
+// T0Scramble() //
+//////////////////
+
+bool T0MTHistos :: T0Scramble()
+ {
+ std::cout<<"Scrambling for "<<m_time->GetName()<<std::endl;
+ string fitopt("BLR");
+ if(m_settings->VerboseLevel()==0)
+ fitopt+="Q";
+ if(m_settings->AddFitfun())
+ {
+ fitopt+="+";
+ }
+ else
+ {
+ fitopt+="0";
+ }
+//create scrambled histogram
+ char scramhistname[100];
+ snprintf(scramhistname, 100, "%s_scram", m_time->GetName());
+ TH1F *scramhist = new TH1F(scramhistname, "scrambled histogram", m_time->GetSize()-2, m_time->GetXaxis()->GetXmin(), m_time->GetXaxis()->GetXmax());
+ for (int binnr=0; binnr < m_time->GetSize(); binnr++)
+ {
+ scramhist->SetBinContent(binnr, m_time->GetBinContent(binnr) + gRandom->Gaus(0,m_time->GetBinError(binnr)));
+ scramhist->SetBinError(binnr, m_time->GetBinError(binnr)*1.41421356);
+ if (scramhist->GetBinContent(binnr)<0) scramhist->SetBinContent(binnr, 0);
+ }
+ TDirectory *cwd=gDirectory;
+ if(dir!=NULL) dir->cd();
+ MTT0PatternRecognition scramrec;
+//perform pattern recognition
+ if(!scramrec.Initialize(scramhist, m_settings))
+ {
+ cerr<<"T0MTHistos :: FitT0 for tube "<<m_id<<": Scrambed pattern recognition failed!"<<endl;
+ cwd->cd();
+ return false;
+ }
+ char scrambuffer[100];
+ snprintf(scrambuffer,100, "scrammt_t0_fermi");
+ TF1 *scramm_t0_fermi=new TF1(scrambuffer, m_t0_fermi, scramrec.GetFitRangeMin(), scramrec.GetFitRangeMax(), N_T0_FIT_PAR);
+//set parameter names
+ scramm_t0_fermi->SetParName(T0_PAR_NR_T0, "t_{0}");
+ scramm_t0_fermi->SetParName(T0_PAR_NR_T, "T");
+ scramm_t0_fermi->SetParName(T0_PAR_NR_BACK, "r_{b}");
+ scramm_t0_fermi->SetParName(T0_PAR_NR_A, "A");
+//set fixed values
+ scramm_t0_fermi->FixParameter(T0_PAR_NR_BACK, scramrec.GetBackground());
+ scramm_t0_fermi->SetParameter(T0_PAR_NR_A, scramrec.GetHeight() - scramrec.GetBackground());
+//set estimates as start values
+ scramm_t0_fermi->SetParameter(T0_PAR_NR_T0, scramrec.GetEstimatedT0());
+//set resonable start value for T
+ scramm_t0_fermi->SetParameter(T0_PAR_NR_T, 3.0);
+//perform fit - NOTE: The return value of the Fit function is not documented!
+ scramhist->Fit(scrambuffer, fitopt.c_str(),"", scramrec.GetFitRangeMin(), scramrec.GetFitRangeMax());
+//set parameter for the new fit of the original histogram
+ m_time->GetListOfFunctions()->Clear();
+//set fixed values
+ m_t0_fermi->FixParameter(T0_PAR_NR_BACK, scramm_t0_fermi->GetParameter(T0_PAR_NR_BACK));
+ m_t0_fermi->SetParameter(T0_PAR_NR_A, scramm_t0_fermi->GetParameter(T0_PAR_NR_A));
+//set estimates as start values
+ m_t0_fermi->SetParameter(T0_PAR_NR_T0, scramm_t0_fermi->GetParameter(T0_PAR_NR_T0));
+//set resonable start value for T
+ m_t0_fermi->SetParameter(T0_PAR_NR_T, scramm_t0_fermi->GetParameter(T0_PAR_NR_T));
+ m_time->GetListOfFunctions()->Clear();
+ m_time->Fit("mt_t0_fermi", fitopt.c_str(),"", scramrec.GetFitRangeMin(), scramrec.GetFitRangeMax());
+ if(m_settings->T0Settings()->UseTopChi2())
+ TopChi2();
+ else
+ m_chi2=m_time->GetFunction("mt_t0_fermi")->GetChisquare() / m_time->GetFunction("mt_t0_fermi")->GetNDF();
+ cwd->cd();
+ return true;
+ }
+
+void T0MTHistos :: TopChi2()
+ {
+//calculate topchi2
+ m_chi2 = 0;
+ int topndf = 0;
+ TF1 *t0_fermi=m_time->GetFunction("mt_t0_fermi");
+ Double_t min, max;
+ t0_fermi->GetRange(min, max);
+ int startbin = m_time->FindBin(min) + 1;
+ int endbin = m_time->FindBin(max) -1;
+ for (int bin = startbin; bin < endbin; bin++)
+ {
+ float measval = m_time->GetBinContent(bin);
+ float funcval = t0_fermi->Eval(m_time->GetBinCenter(bin));
+ float errval = m_time->GetBinError(bin);
+ //take only chi2 from top part or if the bin content is min 10 or if the function >10
+ if(measval < 10 && funcval<10 && m_time->GetBinCenter(bin) < t0_fermi->GetParameter(T0_PAR_NR_T0) + 2 * t0_fermi->GetParameter(T0_PAR_NR_T))
+ continue;
+ if (errval == 0) m_chi2+= ( measval - funcval )*(measval - funcval);
+ else m_chi2+= ( measval - funcval )*(measval - funcval)/(errval*errval);
+ topndf++;
+ }
+ m_chi2=m_chi2/topndf;
+ }
+
+
+void T0MTHistos :: TopSlicing()
+ {
+ std::cout<<"Slicing for "<<m_time->GetName()<<std::endl;
+//vector with slice chi2
+ std::list<Slice> slice_chi2;
+ Slice current;
+ current.chi_2=0.0;
+ current.n_bins=0;
+ TF1 *t0_fermi=m_time->GetFunction("mt_t0_fermi");
+ current.min_bin=m_time->FindBin(t0_fermi->GetParameter(T0_PAR_NR_T0) + 2 * t0_fermi->GetParameter(T0_PAR_NR_T));
+ Double_t min, max;
+ t0_fermi->GetRange(min, max);
+ std::cout<<current.min_bin<<" "<<max<<std::endl;
+ for(int bin=m_time->FindBin(t0_fermi->GetParameter(T0_PAR_NR_T0) + 2 * t0_fermi->GetParameter(T0_PAR_NR_T)); bin<m_time->FindBin(max) - 1; bin++)
+ {
+ if(current.n_bins==10)
+ {
+ std::cout<<current.chi_2/current.n_bins<<std::endl;
+ current.max_bin=bin;
+ slice_chi2.push_back(current);
+ current.chi_2=0.0;
+ current.min_bin=bin;
+ current.n_bins=0;
+ }
+ double measval = m_time->GetBinContent(bin);
+ double funcval = t0_fermi->Eval(m_time->GetBinCenter(bin));
+ double errval = m_time->GetBinError(bin);
+ if(errval==0)
+ current.chi_2 += std::pow(measval - funcval, 2.0);
+ else
+ current.chi_2 += std::pow(measval - funcval, 2.0)/std::pow(errval, 2);
+ current.n_bins++;
+ }
+ std::cout<<"number of slices: "<<slice_chi2.size()<<std::endl;
+ std::list<Slice> :: iterator it=slice_chi2.end();
+ do
+ {
+ it--;
+ if(it==slice_chi2.begin())
+ {
+ std::cout<<"No gain in slicing!"<<std::endl;
+ return;
+ }
+ }
+ while(it->chi_2/static_cast<double>(it->n_bins) > 3);
+ max=m_time->GetBinCenter(it->min_bin);
+ m_time->GetListOfFunctions()->Clear();
+ string fitopt("BLR");
+ if(m_settings->VerboseLevel()==0)
+ fitopt+="Q";
+ if(m_settings->AddFitfun())
+ {
+ fitopt+="+";
+ }
+ else
+ {
+ fitopt+="0";
+ }
+ m_time->Fit("mt_t0_fermi", fitopt.c_str(),"", min, max);
+ if(m_settings->T0Settings()->UseTopChi2())
+ TopChi2();
+ else
+ m_chi2=m_time->GetFunction("mt_t0_fermi")->GetChisquare() / m_time->GetFunction("mt_t0_fermi")->GetNDF();
+ }
+
+
+}
diff --git a/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibT0/src/T0MTSettings.cxx b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibT0/src/T0MTSettings.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..4d92d2e89b4add3f2971c40fcc2ff85bd5b87e93
--- /dev/null
+++ b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibT0/src/T0MTSettings.cxx
@@ -0,0 +1,28 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+
+#include "MdtCalibT0/T0MTSettings.h"
+
+namespace MuonCalib {
+
+T0MTSettings :: T0MTSettings(): m_n_bins_adc(1000),
+ m_adc_min(-0.5),
+ m_adc_max(999.5),
+ m_n_bins_time(1722),
+ m_time_min(-15.625),
+ m_time_max(1330.08),
+ m_draw_debug_graphs(false),
+ m_add_fitfun(false),
+ m_verbose_level(0),
+ m_fit_time(true),
+ m_fit_adc(false),
+ m_min_entries_time(10000),
+ m_min_entries_adc(1000),
+ m_t0_settings(),
+ m_tmax_settings()
+
+ {
+ }
+}
diff --git a/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibT0/src/T0MTSettingsT0.cxx b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibT0/src/T0MTSettingsT0.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..8faa2d99796a88a02cad1e80e543ae5a1e96a9c6
--- /dev/null
+++ b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibT0/src/T0MTSettingsT0.cxx
@@ -0,0 +1,20 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+
+#include "MdtCalibT0/T0MTSettingsT0.h"
+
+namespace MuonCalib {
+
+T0MTSettingsT0 :: T0MTSettingsT0() : m_vbh_bin_content_rel(2.0),
+ m_max_bin_width(10.0),
+ m_min_background_bins(50),
+ m_T_start(3.0),
+ m_use_top_chi2(true),
+ m_scramble_threshold(-1.0),
+ m_slicing_chi2(-1.0),
+ m_correct_rel_t0s(false)
+ {
+ }
+}
diff --git a/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibT0/src/T0MTSettingsTMax.cxx b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibT0/src/T0MTSettingsTMax.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..331e0e39d07fb3231fbb6b3ca2080d3ea65411b3
--- /dev/null
+++ b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibT0/src/T0MTSettingsTMax.cxx
@@ -0,0 +1,23 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+
+#include "MdtCalibT0/T0MTSettingsTMax.h"
+
+namespace MuonCalib {
+
+T0MTSettingsTMax :: T0MTSettingsTMax() : m_vbh_bin_content_rel(4.0),
+ m_max_bin_width(40.0),
+ m_vbh_low(200.0),
+ m_end_min(500.0),
+ m_end_max(1000.0),
+ m_dist_background(80.0),
+ m_min_background_bins(10),
+ m_dist_ab(50.0),
+ m_width_ab(300.0),
+ m_start_t(3.0)
+ {
+ }
+
+}
diff --git a/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibT0/src/VariableBinwidthHistogram.cxx b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibT0/src/VariableBinwidthHistogram.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..dee319077a153daec65ef42cd0054761db34be3c
--- /dev/null
+++ b/MuonSpectrometer/MuonCalib/MdtCalib/MdtCalibT0/src/VariableBinwidthHistogram.cxx
@@ -0,0 +1,220 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+using namespace std;
+
+#include "MdtCalibT0/VariableBinwidthHistogram.h"
+
+// ROOT
+#include "TH1.h"
+#include "TGraph.h"
+#include <cmath>
+
+namespace MuonCalib{
+
+
+//////////////////////////////////////////////////////////////////////////
+// Initialize(const &TH1 hist, double binc_rate, double max_bin_width) //
+//////////////////////////////////////////////////////////////////////////
+
+bool VariableBinwidthHistogram :: Initialize(TH1F * hist, double binc_rate, double max_bin_width, double min_x, double max_x)
+ {
+ if(binc_rate<1.0)
+ {
+ cerr<<"ERROR in VariableBinwidthHistogram :: Initialize!"<<endl;
+ cerr<<"binc_rate must be greater than 1!"<<endl;
+ m_error=true;
+ return false;
+ }
+//get maximum, firest bin and last bin
+ double max;
+ int first_bin(hist->FindBin(min_x)), last_bin(hist->FindBin(max_x));
+ bool max_valid(true);
+ if(first_bin==0) first_bin=1;
+ else max_valid=false;
+ if(last_bin>hist->GetNbinsX())last_bin=hist->GetNbinsX();
+ else max_valid=false;
+ if(max_valid) max=hist->GetMaximum();
+ else
+ {
+ max=-9e9;
+ for(int i=first_bin; i<=last_bin; i++)
+ {
+ if(max<hist->GetBinContent(i)) max=hist->GetBinContent(i);
+ }
+ }
+//get number of entries per bin
+ m_binc=static_cast<unsigned int>(ceil(max*binc_rate));
+ m_max_bin_width=max_bin_width;
+//create first bin
+ m_bins.clear();
+ m_bins.push_back(new VariableBinwidthHistogramBin(hist->GetBinLowEdge(first_bin) + 0.5 * max_bin_width, max_bin_width, 0));
+ VariableBinwidthHistogramBin *current_bin(m_bins[0]);
+//loop on histogram bins
+ for(int i=first_bin; i<=last_bin; i++)
+ {
+ //maximum binwidth reached
+ if(hist->GetBinCenter(i) > current_bin->Right())
+ {
+ m_bins. push_back(new VariableBinwidthHistogramBin(current_bin->Right() + 0.5 * m_max_bin_width, max_bin_width, 0));
+ current_bin=m_bins[m_bins.size()-1];
+ }
+ //will the bin be full
+ if(current_bin->Entries() + static_cast<unsigned int>(hist->GetBinContent(i)) > m_binc)
+ {
+ //this will be filled in the next bin
+ unsigned int overflow=current_bin->Entries() + static_cast<unsigned int>(hist->GetBinContent(i)) - m_binc;
+ //the current bin ends here.
+ current_bin->MoveRight(hist->GetBinLowEdge(i));
+ current_bin->SetContent(m_binc);
+ //create new bin
+ m_bins. push_back(new VariableBinwidthHistogramBin(current_bin->Right() + 0.5 * m_max_bin_width, m_max_bin_width, overflow));
+ current_bin=m_bins[m_bins.size()-1];
+ continue;
+ }
+ //add to content of current bin
+ (*current_bin)+=static_cast<unsigned int>(hist->GetBinContent(i));
+ }
+//sort bins by content
+ for(unsigned int i=0; i<m_bins.size(); i++)
+ {
+ m_sort_bins.push_back(VBHBinPtrSrt(m_bins[i]));
+ }
+ m_sorted=false;
+ m_error=false;
+ return true;
+ }
+
+
+//////////////////////////
+// Smooth(double width) //
+//////////////////////////
+
+bool VariableBinwidthHistogram :: Smooth(double width)
+ {
+//needs at last 3 bins to smooth
+ if(m_bins.size()<3)
+ {
+ cerr<<"VariableBinwidthHistogram :: Smooth: VBH has less than 3 bins!"<<endl;
+ return false;
+ }
+ for(unsigned int i=0; i<m_bins.size()-3; i++)
+ {
+ Double_t sl1=m_bins[i+1]->Width()-m_bins[i]->Width();
+ Double_t sl2=m_bins[i+2]->Width()-m_bins[i+1]->Width();
+ //one slopes must be smaller or equal to bw
+// if(fabs(sl1)>width && fabs(sl2)>width) continue;
+ //slopes must be oposit sign
+ if(sign(sl1)==sign(sl2)) continue;
+ //prevents numerical effects
+ if(fabs(sl1)<width/2 || fabs(sl2)<width/2) continue;
+ //move bin boarder
+ m_bins[i]->MoveRight(m_bins[i]->Right()-width/2*sign(sl2));
+ m_bins[i+1]->MoveLeft(m_bins[i]->Right()-width/2*sign(sl2));
+ }
+ m_sorted=false;
+ return true;
+ }
+
+//////////////////////////
+// DenistyGraph() //
+//////////////////////////
+
+TGraph *VariableBinwidthHistogram :: DenistyGraph() const
+ {
+ Double_t *x = new Double_t[m_bins.size()],
+ *y = new Double_t[m_bins.size()];
+ for(unsigned int i=0; i<m_bins.size(); i++)
+ {
+ x[i]=m_bins[i]->Center();
+ y[i]=m_bins[i]->Density();
+ }
+ TGraph *gr=new TGraph(m_bins.size(), x, y);
+// cout<<"delete"<<endl;
+// delete [] x;
+// delete [] y;
+// cout<<"return"<<endl;
+ return gr;
+ }
+
+
+//////////////////////////////////
+// BinWidthGraph() const //
+//////////////////////////////////
+
+TGraph *VariableBinwidthHistogram :: BinWidthGraph() const
+ {
+ Double_t *x = new Double_t[m_bins.size()],
+ *y = new Double_t[m_bins.size()];
+ for(unsigned int i=0; i<m_bins.size(); i++)
+ {
+ x[i]=m_bins[i]->Center();
+ y[i]=m_bins[i]->Width();
+ }
+ TGraph *gr=new TGraph(m_bins.size(), x, y);
+// delete [] x;
+// delete [] y;
+ return gr;
+ }
+
+
+//////////////////////////////////
+// BinContentGraph() const //
+//////////////////////////////////
+
+TGraph *VariableBinwidthHistogram :: BinContentGraph() const
+ {
+ Double_t *x = new Double_t[m_bins.size()],
+ *y = new Double_t[m_bins.size()];
+ for(unsigned int i=0; i<m_bins.size(); i++)
+ {
+ x[i]=m_bins[i]->Center();
+ y[i]=m_bins[i]->Entries();
+ }
+ TGraph *gr=new TGraph(m_bins.size(), x, y);
+// delete [] x;
+// delete [] y;
+ return gr;
+ }
+
+
+TGraph *VariableBinwidthHistogram :: DiffDensityGraph() const
+ {
+ if(m_bins.size()<2)
+ {
+ cerr<<"VariableBinwidthHistogram :: DiffDensity(): Need at alst 2 bins for differential density!"<<endl;
+ return new TGraph();
+ }
+ Double_t *x = new Double_t[m_bins.size()-1],
+ *y = new Double_t[m_bins.size()-1];
+ for(unsigned int i=0; i<m_bins.size()-1; i++)
+ {
+ x[i]=m_bins[i]->Center();
+ y[i]=m_bins[i+1]->Density() - m_bins[i]->Density();
+ }
+ TGraph *gr=new TGraph(m_bins.size()-1, x, y);
+ return gr;
+ }
+
+
+
+TGraph *VariableBinwidthHistogram :: DiffBinwidthGraph() const
+ {
+ if(m_bins.size()<2)
+ {
+ cerr<<"VariableBinwidthHistogram :: DiffBinwidth(): Need at alst 2 bins for differential density!"<<endl;
+ return new TGraph();
+ }
+ Double_t *x = new Double_t[m_bins.size()-1],
+ *y = new Double_t[m_bins.size()-1];
+ for(unsigned int i=0; i<m_bins.size()-1; i++)
+ {
+ x[i]=m_bins[i]->Center();
+ y[i]=m_bins[i+1]->Width() - m_bins[i]->Width();
+ }
+ TGraph *gr=new TGraph(m_bins.size()-1, x, y);
+ return gr;
+ }
+
+}