MuonResidualsAngleFitter Class Reference

#include <MuonResidualsAngleFitter.h>

Inheritance diagram for MuonResidualsAngleFitter:

Public Types
enum	{   kAngle = 0, kXControl, kYControl, kSigma,   kGamma, kNPar }
enum	{ kResidual = 0, kXAngle, kYAngle, kNData }
Public Types inherited from MuonResidualsFitter
enum	{   kPureGaussian, kPowerLawTails, kROOTVoigt, kGaussPowerTails,   kPureGaussian2D }
enum	{   k1DOF, k5DOF, k6DOF, k6DOFrphi,   kPositionFitter, kAngleFitter, kAngleBfieldFitter }
enum	{   k1111, k1110, k1100, k1010,   k0010, k1000, k0100 }

Public Member Functions
bool	fit (Alignable *ali) override
	MuonResidualsAngleFitter (int residualsModel, int minHitsPerRegion, int useResiduals, bool weightAlignment=true)
int	ndata () override
int	npar () override
double	plot (std::string name, TFileDirectory *dir, Alignable *ali) override
double	sumofweights () override
int	type () const override
Public Member Functions inherited from MuonResidualsFitter
void	computeHistogramRangeAndBinning (int which, int &nbins, double &a, double &b)
virtual void	correctBField ()=0
virtual void	correctBField (int idx_momentum, int idx_q)
TMatrixDSym	correlationMatrix ()
double	covarianceElement (int parNum1, int parNum2)
TMatrixDSym	covarianceMatrix ()
void	eraseNotSelectedResiduals ()
double	errorerror (int parNum)
void	fiducialCuts (double xMin=-80.0, double xMax=80.0, double yMin=-80.0, double yMax=80.0, bool fidcut1=false)
void	fill (double *residual)
void	fix (int parNum, bool dofix=true)
bool	fixed (int parNum)
void	histogramChi2GaussianFit (int which, double &fit_mean, double &fit_sigma)
double	loglikelihood ()
	MuonResidualsFitter (int residualsModel, int minHits, int useResiduals, bool weightAlignment=true)
int	nfixed ()
long	numResiduals () const
long	numsegments ()
int	parNum2parIdx (int parNum)
void	plotsimple (std::string name, TFileDirectory *dir, int which, double multiplier)
void	plotweighted (std::string name, TFileDirectory *dir, int which, int whichredchi2, double multiplier)
void	read (FILE *file, int which=0)
std::vector< double * > ::const_iterator	residuals_begin () const
std::vector< double * > ::const_iterator	residuals_end () const
int	residualsModel () const
std::vector< bool > &	selectedResidualsFlags ()
void	selectPeakResiduals (double nsigma, int nvar, int *vars)
void	selectPeakResiduals_simple (double nsigma, int nvar, int *vars)
void	setInitialValue (int parNum, double value)
void	setPrintLevel (int printLevel)
void	setStrategy (int strategy)
int	useRes (int pattern=-1)
double	value (int parNum)
void	write (FILE *file, int which=0)
virtual	~MuonResidualsFitter ()

Protected Member Functions
void	inform (TMinuit *tMinuit) override
Protected Member Functions inherited from MuonResidualsFitter
bool	dofit (void(*fcn)(int &, double *, double &, double *, int), std::vector< int > &parNum, std::vector< std::string > &parName, std::vector< double > &start, std::vector< double > &step, std::vector< double > &low, std::vector< double > &high)
void	initialize_table ()

Additional Inherited Members
Protected Attributes inherited from MuonResidualsFitter
double	m_center [20]
TMatrixDSym	m_cov
std::vector< double >	m_error
std::vector< bool >	m_fixed
double	m_loglikelihood
int	m_minHits
std::map< int, double >	m_parNum2InitValue
std::map< int, int >	m_parNum2parIdx
int	m_printLevel
double	m_radii [20]
std::vector< double * >	m_residuals
std::vector< bool >	m_residuals_ok
int	m_residualsModel
int	m_strategy
int	m_useResiduals
std::vector< double >	m_value
bool	m_weightAlignment

Detailed Description

Date:

2010/03/12 22:23:18

Revision:

1.10

Author

J. Pivarski - Texas A&M University pivar.nosp@m.ski@.nosp@m.physi.nosp@m.cs.t.nosp@m.amu.e.nosp@m.du

Definition at line 12 of file MuonResidualsAngleFitter.h.

Member Enumeration Documentation

anonymous enum

Enumerator
kAngle
kXControl
kYControl
kSigma
kGamma
kNPar

Definition at line 14 of file MuonResidualsAngleFitter.h.

{ kAngle = 0, kXControl, kYControl, kSigma, kGamma, kNPar };

anonymous enum

Enumerator
kResidual
kXAngle
kYAngle
kNData

Definition at line 16 of file MuonResidualsAngleFitter.h.

{ kResidual = 0, kXAngle, kYAngle, kNData };

Constructor & Destructor Documentation

MuonResidualsAngleFitter::MuonResidualsAngleFitter ( int  residualsModel, int  minHitsPerRegion, int  useResiduals, bool  weightAlignment = true  )

inline

Definition at line 18 of file MuonResidualsAngleFitter.h.

      : MuonResidualsFitter(residualsModel, minHitsPerRegion, useResiduals, weightAlignment) {}

Member Function Documentation

bool MuonResidualsAngleFitter::fit ( Alignable *  ali )

overridevirtual

Implements MuonResidualsFitter.

Definition at line 40 of file MuonResidualsAngleFitter.cc.

References MuonResidualsFitter::dofit(), MuonResidualsFitter::initialize_table(), kAngle, kGamma, MuonResidualsFitter::kGaussPowerTails, MuonResidualsFitter::kPureGaussian, kResidual, kSigma, kXControl, kYControl, MuonResidualsFitter::m_minHits, SiStripPI::mean, MuonResidualsAngleFitter_FCN(), N, funct::pow(), MuonResidualsFitter::residuals_begin(), MuonResidualsFitter::residuals_end(), MuonResidualsFitter::residualsModel(), mathSSE::sqrt(), command_line::start, plotscripts::stdev(), AlCaHLTBitMon_QueryRunRegistry::string, makeMuonMisalignmentScenario::sum_x, and makeMuonMisalignmentScenario::sum_xx.

Referenced by trackingPlots.Iteration::modules().

                                                 {
  initialize_table();  // if not already initialized

  double sum_x = 0.;
  double sum_xx = 0.;
  int N = 0;

  for (std::vector<double *>::const_iterator resiter = residuals_begin(); resiter != residuals_end(); ++resiter) {
    const double residual = (*resiter)[kResidual];
    //     const double xangle = (*resiter)[kXAngle];
    //     const double yangle = (*resiter)[kYAngle];

    if (fabs(residual) < 0.1) {  // truncate at 100 mrad
      sum_x += residual;
      sum_xx += residual * residual;
      N++;
    }
  }

  if (N < m_minHits)
    return false;

  // truncated mean and stdev to seed the fit
  double mean = sum_x / double(N);
  double stdev = sqrt(sum_xx / double(N) - pow(sum_x / double(N), 2));

  // refine the standard deviation calculation
  sum_x = 0.;
  sum_xx = 0.;
  N = 0;
  for (std::vector<double *>::const_iterator resiter = residuals_begin(); resiter != residuals_end(); ++resiter) {
    const double residual = (*resiter)[kResidual];
    if (mean - 1.5 * stdev < residual && residual < mean + 1.5 * stdev) {
      sum_x += residual;
      sum_xx += residual * residual;
      N++;
    }
  }
  mean = sum_x / double(N);
  stdev = sqrt(sum_xx / double(N) - pow(sum_x / double(N), 2));

  sum_x = 0.;
  sum_xx = 0.;
  N = 0;
  for (std::vector<double *>::const_iterator resiter = residuals_begin(); resiter != residuals_end(); ++resiter) {
    const double residual = (*resiter)[kResidual];
    if (mean - 1.5 * stdev < residual && residual < mean + 1.5 * stdev) {
      sum_x += residual;
      sum_xx += residual * residual;
      N++;
    }
  }
  mean = sum_x / double(N);
  stdev = sqrt(sum_xx / double(N) - pow(sum_x / double(N), 2));

  std::vector<int> parNum;
  std::vector<std::string> parName;
  std::vector<double> start;
  std::vector<double> step;
  std::vector<double> low;
  std::vector<double> high;

  parNum.push_back(kAngle);
  parName.push_back(std::string("angle"));
  start.push_back(mean);
  step.push_back(0.1);
  low.push_back(0.);
  high.push_back(0.);
  parNum.push_back(kXControl);
  parName.push_back(std::string("xcontrol"));
  start.push_back(0.);
  step.push_back(0.1);
  low.push_back(0.);
  high.push_back(0.);
  parNum.push_back(kYControl);
  parName.push_back(std::string("ycontrol"));
  start.push_back(0.);
  step.push_back(0.1);
  low.push_back(0.);
  high.push_back(0.);
  parNum.push_back(kSigma);
  parName.push_back(std::string("sigma"));
  start.push_back(stdev);
  step.push_back(0.1 * stdev);
  low.push_back(0.);
  high.push_back(0.);
  if (residualsModel() != kPureGaussian && residualsModel() != kGaussPowerTails) {
    parNum.push_back(kGamma);
    parName.push_back(std::string("gamma"));
    start.push_back(stdev);
    step.push_back(0.1 * stdev);
    low.push_back(0.);
    high.push_back(0.);
  }

  return dofit(&MuonResidualsAngleFitter_FCN, parNum, parName, start, step, low, high);
}

void MuonResidualsAngleFitter::inform ( TMinuit *  tMinuit )

overrideprotectedvirtual

Implements MuonResidualsFitter.

Definition at line 5 of file MuonResidualsAngleFitter.cc.

References MuonResidualsAngleFitter_TMinuit.

{ MuonResidualsAngleFitter_TMinuit = tMinuit; }

int MuonResidualsAngleFitter::ndata ( )

inlineoverridevirtual

Implements MuonResidualsFitter.

Definition at line 33 of file MuonResidualsAngleFitter.h.

References kNData.

{ return kNData; }

int MuonResidualsAngleFitter::npar ( )

inlineoverridevirtual

Implements MuonResidualsFitter.

Definition at line 23 of file MuonResidualsAngleFitter.h.

References cms::cuda::assert(), MuonResidualsFitter::kGaussPowerTails, kNPar, MuonResidualsFitter::kPowerLawTails, MuonResidualsFitter::kPureGaussian, MuonResidualsFitter::kROOTVoigt, and MuonResidualsFitter::residualsModel().

                      {
    if (residualsModel() == kPureGaussian || residualsModel() == kGaussPowerTails)
      return kNPar - 1;
    else if (residualsModel() == kPowerLawTails)
      return kNPar;
    else if (residualsModel() == kROOTVoigt)
      return kNPar;
    else
      assert(false);
  }

double MuonResidualsAngleFitter::plot ( std::string  name, TFileDirectory *  dir, Alignable *  ali  )

overridevirtual

Implements MuonResidualsFitter.

Definition at line 138 of file MuonResidualsAngleFitter.cc.

References kAngle, kGamma, MuonResidualsFitter::kGaussPowerTails, MuonResidualsFitter::kPowerLawTails, MuonResidualsFitter::kPureGaussian, kResidual, MuonResidualsFitter::kROOTVoigt, kSigma, kXAngle, kXControl, kYAngle, kYControl, TFileDirectory::make(), MuonResidualsFitter_GaussPowerTails_TF1(), MuonResidualsFitter_powerLawTails_TF1(), MuonResidualsFitter_pureGaussian_TF1(), MuonResidualsFitter_ROOTVoigt_TF1(), MuonResidualsFitter::numResiduals(), MuonResidualsFitter::residuals_begin(), MuonResidualsFitter::residuals_end(), MuonResidualsFitter::residualsModel(), AlCaHLTBitMon_QueryRunRegistry::string, and MuonResidualsFitter::value().

                                                                                         {
  std::stringstream raw_name, narrowed_name, xcontrol_name, ycontrol_name;
  raw_name << name << "_raw";
  narrowed_name << name << "_narrowed";
  xcontrol_name << name << "_xcontrol";
  ycontrol_name << name << "_ycontrol";

  TH1F *raw_hist =
      dir->make<TH1F>(raw_name.str().c_str(), (raw_name.str() + std::string(" (mrad)")).c_str(), 100, -100., 100.);
  TH1F *narrowed_hist = dir->make<TH1F>(
      narrowed_name.str().c_str(), (narrowed_name.str() + std::string(" (mrad)")).c_str(), 100, -100., 100.);
  TProfile *xcontrol_hist = dir->make<TProfile>(
      xcontrol_name.str().c_str(), (xcontrol_name.str() + std::string(" (mrad)")).c_str(), 100, -1., 1.);
  TProfile *ycontrol_hist = dir->make<TProfile>(
      ycontrol_name.str().c_str(), (ycontrol_name.str() + std::string(" (mrad)")).c_str(), 100, -1., 1.);

  narrowed_name << "fit";
  xcontrol_name << "fit";
  ycontrol_name << "fit";

  double scale_factor = double(numResiduals()) * (100. - -100.) / 100;  // (max - min)/nbins

  TF1 *narrowed_fit = nullptr;
  if (residualsModel() == kPureGaussian) {
    narrowed_fit = new TF1(narrowed_name.str().c_str(), MuonResidualsFitter_pureGaussian_TF1, -100., 100., 3);
    narrowed_fit->SetParameters(scale_factor, value(kAngle) * 1000., value(kSigma) * 1000.);
    narrowed_fit->Write();
  } else if (residualsModel() == kPowerLawTails) {
    narrowed_fit = new TF1(narrowed_name.str().c_str(), MuonResidualsFitter_powerLawTails_TF1, -100., 100., 4);
    narrowed_fit->SetParameters(scale_factor, value(kAngle) * 1000., value(kSigma) * 1000., value(kGamma) * 1000.);
    narrowed_fit->Write();
  } else if (residualsModel() == kROOTVoigt) {
    narrowed_fit = new TF1(narrowed_name.str().c_str(), MuonResidualsFitter_ROOTVoigt_TF1, -100., 100., 4);
    narrowed_fit->SetParameters(scale_factor, value(kAngle) * 1000., value(kSigma) * 1000., value(kGamma) * 1000.);
    narrowed_fit->Write();
  } else if (residualsModel() == kGaussPowerTails) {
    narrowed_fit = new TF1(narrowed_name.str().c_str(), MuonResidualsFitter_GaussPowerTails_TF1, -100., 100., 3);
    narrowed_fit->SetParameters(scale_factor, value(kAngle) * 1000., value(kSigma) * 1000.);
    narrowed_fit->Write();
  }

  TF1 *xcontrol_fit = new TF1(xcontrol_name.str().c_str(), "[0]+x*[1]", -1., 1.);
  xcontrol_fit->SetParameters(value(kAngle) * 1000., value(kXControl) * 1000.);
  xcontrol_fit->Write();

  TF1 *ycontrol_fit = new TF1(ycontrol_name.str().c_str(), "[0]+x*[1]", -1., 1.);
  ycontrol_fit->SetParameters(value(kAngle) * 1000., value(kYControl) * 1000.);
  ycontrol_fit->Write();

  for (std::vector<double *>::const_iterator resiter = residuals_begin(); resiter != residuals_end(); ++resiter) {
    const double raw_residual = (*resiter)[kResidual];
    const double xangle = (*resiter)[kXAngle];
    const double yangle = (*resiter)[kYAngle];

    double xangle_correction = value(kXControl) * xangle;
    double yangle_correction = value(kYControl) * yangle;
    double corrected_residual = raw_residual - xangle_correction - yangle_correction;

    raw_hist->Fill(raw_residual * 1000.);
    narrowed_hist->Fill(corrected_residual * 1000.);

    xcontrol_hist->Fill(xangle, (raw_residual - yangle_correction) * 1000.);
    ycontrol_hist->Fill(yangle, (raw_residual - xangle_correction) * 1000.);
  }

  return 0.;
}

double MuonResidualsAngleFitter::sumofweights ( )

inlineoverridevirtual

Implements MuonResidualsFitter.

Definition at line 36 of file MuonResidualsAngleFitter.h.

References MuonResidualsFitter::numResiduals().

{ return numResiduals(); }

int MuonResidualsAngleFitter::type ( ) const

inlineoverridevirtual

Implements MuonResidualsFitter.

Definition at line 21 of file MuonResidualsAngleFitter.h.

References MuonResidualsFitter::kAngleFitter.

{ return MuonResidualsFitter::kAngleFitter; }