MuonCalib::T0CalibrationClassic Class Reference

#include <T0CalibrationClassic.h>

Inheritance diagram for MuonCalib::T0CalibrationClassic:

Collaboration diagram for MuonCalib::T0CalibrationClassic:

Public Types
using	MuonSegVec = std::vector< std::shared_ptr< MuonCalibSegment > >
using	MuonSegIt = MuonSegVec::iterator
using	MuonSegCit = MuonSegVec::const_iterator
using	MdtCalibOutputPtr = std::shared_ptr< IMdtCalibrationOutput >

Public Member Functions
	T0CalibrationClassic (const std::string &name, const T0ClassicSettings *settings)
	constructor More...
	~T0CalibrationClassic ()
	destructor More...
bool	handleSegment (MuonCalibSegment &seg)
	fill tube spectra More...
void	setInput (const IMdtCalibrationOutput *input) override
	unused More...
bool	analyse ()
	extract parameters from spectra More...
bool	converged () const
	return m_converged (always false?) More...
virtual MdtCalibOutputPtr	getResults () const override
virtual MdtCalibOutputPtr	analyseSegments (const MuonSegVec &segs) override
	new interface function More...
virtual std::string	name () const
	returns name (region) of instance More...

Private Member Functions
T0ClassicHistos *	bookHistos (unsigned int idtube)
	booking of histograms More...
T0ClassicHistos *	getHistos (unsigned int idtube)
	retrieve pointer for tube idtube histograms More...
void	doTimeFit (T0ClassicHistos &, MdtTubeFitContainer::SingleTubeFit &, MdtTubeFitContainer::SingleTubeCalib &)
	fit time spectrum More...
void	doAdcFit (T0ClassicHistos &, MdtTubeFitContainer::SingleTubeFit &, MdtTubeFitContainer::SingleTubeCalib &)
	fit adc spectrum More...
void	searchParams (TH1 *h, double *p, int np)
	estimate initial parameters for time spectrum fit from the spectrum itself More...
T0CalibrationClassic &	operator= (const T0CalibrationClassic &right)=delete
	T0CalibrationClassic (const T0CalibrationClassic &)=delete

Private Attributes
const T0ClassicSettings *	m_settings
	pointer to the settings More...
bool	m_converged
	convergence status More...
std::string	m_name
	calibration region name More...
int	m_currentItnum
	current iteration (always 1?) More...
std::unique_ptr< TFile >	m_file {}
	pointer to the histogram file More...
TDirectory *	m_regiondir
	pointer to the ROOT directory More...
std::vector< std::unique_ptr< T0ClassicHistos > >	m_histos
	vector of pointers to tube histograms More...
std::unique_ptr< MdtTubeFitContainer >	m_result
	tube constants More...
bool	m_delete_settings

Detailed Description

Implementation of a T0 calibration using the classical approach.

Definition at line 119 of file T0CalibrationClassic.h.

Member Typedef Documentation

MdtCalibOutputPtr

using MuonCalib::IMdtCalibration::MdtCalibOutputPtr = std::shared_ptr<IMdtCalibrationOutput>

inherited

Definition at line 30 of file IMdtCalibration.h.

MuonSegCit

using MuonCalib::IMdtCalibration::MuonSegCit = MuonSegVec::const_iterator

inherited

Definition at line 29 of file IMdtCalibration.h.

MuonSegIt

using MuonCalib::IMdtCalibration::MuonSegIt = MuonSegVec::iterator

inherited

Definition at line 28 of file IMdtCalibration.h.

MuonSegVec

using MuonCalib::IMdtCalibration::MuonSegVec = std::vector<std::shared_ptr<MuonCalibSegment> >

inherited

Definition at line 27 of file IMdtCalibration.h.

Constructor & Destructor Documentation

T0CalibrationClassic() [1/2]

MuonCalib::T0CalibrationClassic::T0CalibrationClassic	(	const std::string &	name,
		const T0ClassicSettings *	settings
	)

constructor

Parameters

[in]	name	of the region/chamber to be calibrated
[in]	pointer	to settings vector

Definition at line 30 of file T0CalibrationClassic.cxx.

                                                                                                       :
        IMdtCalibration(name), m_settings(settings), m_converged(false), m_name(name), m_result(nullptr), m_delete_settings(false) {
        if (!m_settings) {
            std::array<double, 8> params{};  // 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;
        }
        MsgStream log(Athena::getMessageSvc(), "T0ClassicSettings");
        if (log.level() <= MSG::INFO) log << MSG::INFO << "T0CalibrationClassic::T0CalibrationClassic" << m_name << " " << name << endmsg;
        if (m_settings->printLevel() <= 2) m_settings->print();
        m_currentItnum = 1;
        std::string HistoFileName = "T0Classic_" + m_name + ".root";
        if (log.level() <= MSG::INFO)
            log << MSG::INFO << "T0CalibrationClassic::T0CalibrationClassic" << m_name << " " << name << " " << HistoFileName << endmsg;
        m_file = std::make_unique<TFile>(HistoFileName.c_str(), "recreate");
        m_regiondir = m_file->mkdir(m_name.c_str());
    }

~T0CalibrationClassic()

MuonCalib::T0CalibrationClassic::~T0CalibrationClassic

(

)

destructor

Definition at line 56 of file T0CalibrationClassic.cxx.

                                                {
        m_file->Write();
        m_file->Close();
        if (m_delete_settings) delete m_settings;
    }

T0CalibrationClassic() [2/2]

MuonCalib::T0CalibrationClassic::T0CalibrationClassic ( const T0CalibrationClassic &  )

privatedelete

Member Function Documentation

analyse()

bool MuonCalib::T0CalibrationClassic::analyse

(

)

extract parameters from spectra

Definition at line 119 of file T0CalibrationClassic.cxx.

                                       {
        MsgStream log(Athena::getMessageSvc(), "T0ClassicSettings");
        if (log.level() <= MSG::INFO) log << MSG::INFO << "T0CalibrationClassic::analyse iteration " << m_currentItnum << endmsg;
 
        const MdtIdHelper& idHelper{m_result->idHelperSvc()->mdtIdHelper()};
        // loop over m_histos histograms
        for (std::unique_ptr<T0ClassicHistos> &hist : m_histos) {
            if (m_settings->fitTime()) {
                MdtTubeFitContainer::SingleTubeFit full;
                MdtTubeFitContainer::SingleTubeCalib st;
                int idtube = hist->id;
                if (idtube != 0 && (int)(idtube / 100000000) != 9) {
                    doTimeFit(*hist, full, st);
                    doAdcFit(*hist, full, st);
                    MuonFixedId fId(idtube);
                    int nML = fId.mdtMultilayer();
                    int nL = fId.mdtTubeLayer();
                    int nT = fId.mdtTube();
                    const Identifier tubeId = idHelper.channelID(fId.stationNameString(), fId.eta(), fId.phi(), 
                                                                nML, nL, nT);
               
                    bool setInfo = m_result->setCalib(st, tubeId, log);
                    if (!setInfo) log << MSG::WARNING << "T0CalibrationClassic::PROBLEM! could not set SingleTubeCalib info " << endmsg;
                    setInfo = m_result->setFit(std::move(full), tubeId, log);
                    if (!setInfo) log << MSG::WARNING << "T0CalibrationClassic::PROBLEM! could not set SingleTubeFullInfo info " << endmsg;
                }
            }
        }
 
        m_currentItnum++;
        m_file->Write();
        return true;
    }

analyseSegments()

T0CalibrationClassic::MdtCalibOutputPtr MuonCalib::T0CalibrationClassic::analyseSegments ( const MuonSegVec &  segs )

overridevirtual

new interface function

Implements MuonCalib::IMdtCalibration.

Definition at line 153 of file T0CalibrationClassic.cxx.

                                                                                                      {
        for (const auto & seg : segs) handleSegment(*seg);
        analyse();
        return getResults();
    }

bookHistos()

T0ClassicHistos * MuonCalib::T0CalibrationClassic::bookHistos ( unsigned int  idtube )

private

booking of histograms

Definition at line 581 of file T0CalibrationClassic.cxx.

                                                                         {
        std::unique_ptr<T0ClassicHistos> histos = std::make_unique<T0ClassicHistos>();
        std::string histonametdc{}, histonameadc{};
 
        histos->id = idtube;
        m_regiondir->cd();
        if ((int)(idtube / 100000000) == 9) {
            int mezz = (idtube) % (900000000);
            histonametdc = std::format("time_mezz_{:}", mezz);
            histonameadc = std::format("charge_mezz_{:}", 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 = std::format("time_{:}_{:}_{:}_{:}", stationName, eta, phi, tubeid);
            histonameadc = std::format("charge_{:}_{:}_{:}_{:}", stationName, eta, phi, tubeid);
        }
 
        histos->time =
            std::make_unique<TH1F>(histonametdc.c_str(), "Drift Time", m_settings->binTime(), m_settings->minTime(), m_settings->maxTime());
        histos->adc = std::make_unique<TH1F>(histonameadc.c_str(), "ADC", m_settings->binAdc(), m_settings->minAdc(), m_settings->maxAdc());
 
        m_histos.emplace_back(std::move(histos));
        return m_histos.back().get();
    }

converged()

bool MuonCalib::T0CalibrationClassic::converged

(

)

const

return m_converged (always false?)

Definition at line 654 of file T0CalibrationClassic.cxx.

{ return m_converged; }

doAdcFit()

void MuonCalib::T0CalibrationClassic::doAdcFit ( T0ClassicHistos &  T0h, MdtTubeFitContainer::SingleTubeFit &  fi, MdtTubeFitContainer::SingleTubeCalib &  stc  )

private

fit adc spectrum

Definition at line 327 of file T0CalibrationClassic.cxx.

                                                                                 {
        // 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);
        std::array<double, 4> par{0}, errpar{0};
 
        MsgStream log(Athena::getMessageSvc(), "T0ClassicSettings");
        MuonFixedId fId(T0h.id);
        if (fId.isValid() && fId.is_mdt()) {
            if (log.level() <= MSG::DEBUG) log << MSG::DEBUG << " doAdcFit : checking Single tube entries" << endmsg;
            double adcThreshold = 50.;
            TH1 *hcheck = T0h.adc.get();
            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());
            if (log.level() <= MSG::DEBUG)
                log << MSG::DEBUG << " doAdcFit : TotHits, nhitsAboveAdcCut " << nhits << " " << nhitsAboveAdcCut << endmsg;
 
            fi.cov[20] = nhits;
            fi.cov[21] = nhitsAboveAdcCut;
 
            if (log.level() <= MSG::DEBUG) log << MSG::DEBUG << " STARTING doAdcFit " << endmsg;
 
            TH1 *h = nullptr;
 
            if (fitMezz == 1) {
                int hIdMezz = fId.mdtMezzanine();
                std::string HistoId = std::format("charge_mezz_{:}", (hIdMezz) % (900000000));
                if (log.level() <= MSG::DEBUG) {
                    log << MSG::DEBUG << " doAdcFit HistogramId : " << T0h.id << endmsg;
                    log << MSG::DEBUG << " doAdcFit Histogram : " << HistoId << endmsg;
                }
                TH1F *adcHis = (TH1F *)m_regiondir->Get(HistoId.c_str());
                if (!adcHis) return;
 
                T0ClassicHistos *histosMezz = getHistos(fId.mdtMezzanine());
                h = histosMezz->adc.get();
            }
 
            if (log.level() <= MSG::VERBOSE)
                log << MSG::VERBOSE << " histogram " << h->GetName() << " " << h->GetTitle() << " entries=" << h->GetEntries() << endmsg;
            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();
                    std::array<double, 4> adcpar{};
                    adcpar[0] = maxval * 2.;
                    adcpar[1] = m;
                    adcpar[2] = r;
                    adcpar[3] = r / 3.;
 
                    std::unique_ptr<TF1> AdcSpectrum = std::make_unique<TF1>(
                        "AdcSpectrum", " ([0]*exp((x-[1])/[2]))/ (1.+exp((x-[1])/[3])) ", m_settings->minAdc(), m_settings->maxAdc());
                    AdcSpectrum->SetParameters(adcpar.data());
                    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 (log.level() <= MSG::VERBOSE)
                        log << MSG::VERBOSE << "chi2/ndof=" << chi2 / ndof << " "
                            << "Mean=" << m << " "
                            << "RMS=" << r << " par 0 1 2 3 " << par[0] << " " << par[1] << " " << par[2] << " " << par[3] << endmsg;
                }
                // THE NEW HISTOGRAM HAS BEEN FITTED
            }
 
            stc.adcCal = par[1];
            fi.adc_par[0] = par[1];
            fi.adc_par[1] = par[2];
            fi.adc_err[0] = errpar[1];
            fi.adc_err[1] = errpar[2];
            if (std::abs(ndof) == 0) {
                fi.adc_chi2 = -1;
            } else {
                fi.adc_chi2 = chi2 / static_cast<float>(ndof);
            }
        }
    }

doTimeFit()

void MuonCalib::T0CalibrationClassic::doTimeFit ( T0ClassicHistos &  T0h, MdtTubeFitContainer::SingleTubeFit &  fi, MdtTubeFitContainer::SingleTubeCalib &  stc  )

private

fit time spectrum

Definition at line 159 of file T0CalibrationClassic.cxx.

                                                                                  {
        std::unique_ptr<TMinuit> gMinuit = std::make_unique<TMinuit>();
        int fitMezz(1);
 
        const int np = m_settings->numParams();
        std::vector<double> pfit(np, 0), errfit(np, 0);
        Double_t **matrix = new Double_t *[np];
        // initialize everything to make Coverity happy
        for (int i = 0; i < np; i++) {
            matrix[i] = new Double_t[np];
            for (int j = 0; j < np; j++) matrix[i][j] = 0.;
        }
        const std::array<double, 8> &pdefault = m_settings->params();
        double chi2{0};
        int ndof{0};
 
        MuonFixedId fId(T0h.id);
        int isMultilayer(0);
 
        if (T0h.id == 1 || T0h.id == 2 || (T0h.id > 10 && T0h.id <= 23)) isMultilayer = 1;
 
        MsgStream log(Athena::getMessageSvc(), "T0ClassicSettings");
        if ((fId.isValid() && fId.is_mdt()) || isMultilayer) {
            if (log.level() <= MSG::INFO) log << MSG::INFO << " STARTING doTimeFit " << endmsg;
            TH1 *h = nullptr;
 
            if (isMultilayer) {
                if (log.level() <= MSG::DEBUG) log << MSG::DEBUG << T0h.id << endmsg;
                h = T0h.time.get();
            }
 
            if (fitMezz == 1 && !isMultilayer) {  // FIT MEZZANINE
                int hIdMezz = fId.mdtMezzanine();
 
                std::string HistoId = std::format("time_mezz_{:}", (hIdMezz) % (900000000));
                if (log.level() <= MSG::DEBUG) {
                    log << MSG::DEBUG << " doTimeFit HistogramId : " << T0h.id << endmsg;
                    log << MSG::DEBUG << " doTimeFit Histogram : " << HistoId << endmsg;
                }
                TH1F *timeHis = (TH1F *)m_regiondir->Get(HistoId.c_str());
                if (!timeHis) {
                    for (int i = 0; i < np; i++) delete[] matrix[i];
                    delete[] matrix;
                    return;
                }
 
                T0ClassicHistos *histosMezz = getHistos(fId.mdtMezzanine());
                h = histosMezz->time.get();
            }
 
            Stat_t entries = h->GetEntries();
            fi.statistics = (int)entries;
 
            if (log.level() <= MSG::VERBOSE)
                log << MSG::VERBOSE << " 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
                std::unique_ptr<TF1> FitFunction{h->GetFunction("TimeSpectrum")};
                if (FitFunction) {
                    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
                    std::unique_ptr<TF1> TimeSpectrum = std::make_unique<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());
                    for (int i = 0; i < np; i++) {
                        pfit[i] = pdefault[i];
                        if (log.level() <= MSG::DEBUG)
                            log << MSG::DEBUG << "T0CalibrationClassic::doTimeFit initial parameters" << i << "=" << pfit[i] << endmsg;
                    }
                    searchParams(h, &pfit[0], np);
                    if (log.level() <= MSG::DEBUG) {
                        log << MSG::DEBUG << "T0CalibrationClassic::doTimeFit parameters after searchParams " << endmsg;
                        for (int i = 0; i < np; ++i) { log << MSG::DEBUG << "i,pfit(i) " << i << " " << pfit[i] << endmsg; }
                    }
                    TimeSpectrum->SetParameters(pfit.data());
                    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);
                    }
                    chi2 = TimeSpectrum->GetChisquare();  // total chi2
                    ndof = TimeSpectrum->GetNDF();        // number of degrees of freedom
                }
                // THE NEW HISTOGRAM HAS BEEN FITTED
                if (ndof == 0.) ndof = -1;
 
                if (log.level() <= MSG::VERBOSE)
                    log << MSG::VERBOSE << " fit results chi2/ndof=" << chi2 / ndof << " T0=" << pfit[4] << " err=" << errfit[4] << endmsg;
                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
                /*******************************************************
 
                  *******************************************************/
                // 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
            }
        }
        for (int i = 0; i < np; i++) delete[] matrix[i];
        delete[] matrix;
 
        if (log.level() <= MSG::DEBUG) log << MSG::DEBUG << " ENDING doTimeFit " << endmsg;
    }

getHistos()

T0ClassicHistos * MuonCalib::T0CalibrationClassic::getHistos ( unsigned int  idtube )

private

retrieve pointer for tube idtube histograms

Definition at line 530 of file T0CalibrationClassic.cxx.

                                                                        {
        std::string HistoId;
        if ((int)(idtube / 100000000) == 9) {
            int mezz = (idtube) % (900000000);
            HistoId = std::format("time_mezz_{:}", mezz);
        } else if (idtube == 0) {
            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 {
            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 = std::format("time_{:}_{:}_{:}_{:}", stationName, eta, phi, tubeid);
        }
        T0ClassicHistos *ret = nullptr;
        TH1F *timeHis = (TH1F *)m_regiondir->Get(HistoId.c_str());
        // We will either find the pointer to an existing histogram or create a new one
        if (!timeHis) {  // book histo if it does not exist
            ret = bookHistos(idtube);
        } else {  // else loop over m_histos histograms to look for the set of histos of tube idtube
            for (const std::unique_ptr<T0ClassicHistos> &hist : m_histos) {
                if (hist->time.get() == timeHis) {
                    ret = hist.get();
                    break;
                }
            }
        }
        return ret;
    }  // end T0CalibrationClassic::getHistos

getResults()

IMdtCalibration::MdtCalibOutputPtr MuonCalib::T0CalibrationClassic::getResults ( ) const

overridevirtual

Returns

the calibration results

Implements MuonCalib::IMdtCalibration.

Definition at line 526 of file T0CalibrationClassic.cxx.

                                                                            {
        return std::make_shared<T0CalibrationOutput>(m_result.get());
    }

handleSegment()

bool MuonCalib::T0CalibrationClassic::handleSegment

(

MuonCalibSegment &

seg

)

fill tube spectra

Definition at line 62 of file T0CalibrationClassic.cxx.

                                                                  {
        for (const MuonCalibSegment::MdtHitPtr& hit : seg.mdtHOT()) {
            float distanceToRO = hit->distanceToReadout();
 
            bool ROside = distanceToRO < 130000.;  // this means that there is no selection along the tube
            if (ROside) {
                const MuonFixedId &id = hit->identify();
                const MdtIdHelper& idHelper{m_result->idHelperSvc()->mdtIdHelper()};
                // get the T0 originally subtracted for this hit
                int nML = id.mdtMultilayer();
                int nL = id.mdtTubeLayer();
                int nT = id.mdtTube();
                const Identifier tubeId = idHelper.channelID(id.stationNameString(), id.eta(), id.phi(), 
                                                             nML, nL, nT);
                const MdtTubeFitContainer::SingleTubeCalib *stc = m_result->getCalib(tubeId);
                if (!stc) {
                    MsgStream log(Athena::getMessageSvc(), "T0ClassicSettings");
                    log << MSG::WARNING << "no Single Tube Calib info found for ML=" << nML << " L=" << nL << " T=" << nT << endmsg;
                    log << MSG::WARNING << "container size " << m_result->size() << endmsg;
                    log << MSG::WARNING << "container nML " << m_result->numMultilayers() << endmsg;
                    log << MSG::WARNING << "container nL " << m_result->numLayers() << endmsg;
                    log << MSG::WARNING << "container nT " << m_result->numTubes() << endmsg;
                }
 
                // get histos
                T0ClassicHistos *histos = getHistos(id.getIdInt());
 
                // fill histos
 
                float ttime = hit->driftTime();
 
                histos->time->Fill(ttime);
                histos->adc->Fill(hit->adcCount());
                // book an additional dummy set of histograms to fit the global T0
                T0ClassicHistos *histosAll = getHistos(0);
                histosAll->time->Fill(ttime);
                histosAll->adc->Fill(hit->adcCount());
 
                // M.I. Jun20-07 ---- Adding MultiLayer Histos :
                T0ClassicHistos *histosML = getHistos(nML);
                histosML->time->Fill(ttime);
                histosML->adc->Fill(hit->adcCount());
                // M.I. Jan1507 ---- Adding Mezzanine Histos :
                T0ClassicHistos *histosMezz = getHistos(id.mdtMezzanine());
                histosMezz->time->Fill(ttime);
                histosMezz->adc->Fill(hit->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(hit->adcCount());
            }  // end if (ROside)
        }      // end loop over seg.mdtHOT
        return true;
    }

name()

virtual std::string MuonCalib::IMdtCalibration::name ( ) const

inlinevirtualinherited

returns name (region) of instance

Definition at line 49 of file IMdtCalibration.h.

{ return m_name; }

operator=()

T0CalibrationClassic& MuonCalib::T0CalibrationClassic::operator= ( const T0CalibrationClassic &  right )

privatedelete

searchParams()

void MuonCalib::T0CalibrationClassic::searchParams ( TH1 *  h, double *  p, int  np  )

private

estimate initial parameters for time spectrum fit from the spectrum itself

Definition at line 433 of file T0CalibrationClassic.cxx.

                                                                     {
        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
        std::unique_ptr<TH1> hnew{h->Rebin(RebinFactor, "hnew")};  // creates a new histogram hnew
        // merging 5 bins of h1 in one bin
        MsgStream log(Athena::getMessageSvc(), "T0ClassicSettings");
        if (log.level() <= MSG::DEBUG)
            log << MSG::DEBUG << "nbinsx,sizex,rebinfactor=" << nbinsX << " " << sizeX << " " << RebinFactor << endmsg;
        float minDeriv(9999.);
        int minDerivBin(0);
        sizeX = hnew->GetBinWidth(1);
        int newbins = hnew->GetNbinsX();
        for (int i = 0; i < np; ++i) {
            if (log.level() <= MSG::DEBUG) log << MSG::DEBUG << "i,p(i) " << i << " " << p[i] << endmsg;
        }
        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.; }
        if (log.level() <= MSG::DEBUG) log << MSG::DEBUG << " t0guess is " << t0guess << endmsg;
        //
        // =================== 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);
        if (log.level() <= MSG::DEBUG) log << MSG::DEBUG << " noise is " << noise << endmsg;
        //
        // =================== Normalization =========================================
        //
        int t0bin = minDerivBin;
        int ix1 = t0bin + (int)(50 / sizeX);
        int ix2 = t0bin + (int)(500 / sizeX);
        if (log.level() <= MSG::DEBUG) log << MSG::DEBUG << "t0bin,ix1,ix2 " << t0bin << " " << ix1 << " " << ix2 << endmsg;
        float P1 = p[1];
        float P2 = p[2];
        float P3 = p[3];
        if (log.level() <= MSG::DEBUG) log << MSG::DEBUG << "P1,P2,P3 start are " << P1 << " " << P2 << " " << P3 << endmsg;
        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 = std::exp(-(a1 - t0guess) / P3);
                float A2 = std::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;
                if (log.level() <= MSG::DEBUG) {
                    log << MSG::DEBUG << "a1,a2 " << a1 << " " << a2 << " cont1,cont2 " << cont1 << " " << cont2 << " A1,A2 " << A1 << " "
                        << A2 << endmsg;
                    log << MSG::DEBUG << " t0Guess .... P1, P2 " << P1 << " " << P2 << endmsg;
                }
                p[1] = P1;
                p[2] = P2;
            }
        }
        return;
    }

setInput()

void MuonCalib::T0CalibrationClassic::setInput ( const IMdtCalibrationOutput *  input )

overridevirtual

unused

Implements MuonCalib::IMdtCalibration.

Definition at line 640 of file T0CalibrationClassic.cxx.

                                                                             {
        // 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);
        if (t0Input) {
            m_result = std::make_unique<MdtTubeFitContainer>(*t0Input->t0s());
            m_result->setImplementation("T0CalibrationClassic");
        }
    }

Member Data Documentation

m_converged

bool MuonCalib::T0CalibrationClassic::m_converged

private

convergence status

Definition at line 146 of file T0CalibrationClassic.h.

m_currentItnum

int MuonCalib::T0CalibrationClassic::m_currentItnum

private

current iteration (always 1?)

Definition at line 148 of file T0CalibrationClassic.h.

m_delete_settings

bool MuonCalib::T0CalibrationClassic::m_delete_settings

private

Definition at line 154 of file T0CalibrationClassic.h.

m_file

std::unique_ptr<TFile> MuonCalib::T0CalibrationClassic::m_file {}

private

pointer to the histogram file

Definition at line 149 of file T0CalibrationClassic.h.

m_histos

std::vector<std::unique_ptr<T0ClassicHistos> > MuonCalib::T0CalibrationClassic::m_histos

private

vector of pointers to tube histograms

Definition at line 152 of file T0CalibrationClassic.h.

m_name

std::string MuonCalib::T0CalibrationClassic::m_name

private

calibration region name

Definition at line 147 of file T0CalibrationClassic.h.

m_regiondir

TDirectory* MuonCalib::T0CalibrationClassic::m_regiondir

private

pointer to the ROOT directory

Definition at line 150 of file T0CalibrationClassic.h.

m_result

std::unique_ptr<MdtTubeFitContainer> MuonCalib::T0CalibrationClassic::m_result

private

tube constants

Definition at line 153 of file T0CalibrationClassic.h.

m_settings

const T0ClassicSettings* MuonCalib::T0CalibrationClassic::m_settings

private

pointer to the settings

Definition at line 145 of file T0CalibrationClassic.h.

---

The documentation for this class was generated from the following files:

• T0CalibrationClassic.h
• T0CalibrationClassic.cxx