MuonResidualsFitter Class Reference

#include <MuonResidualsFitter.h>

Inheritance diagram for MuonResidualsFitter:

Classes
struct	MuonAlignmentTreeRow

Public Types
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
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)
virtual bool	fit (Alignable *ali)=0
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)
virtual int	ndata ()=0
int	nfixed ()
virtual int	npar ()=0
long	numResiduals () const
long	numsegments ()
int	parNum2parIdx (int parNum)
virtual double	plot (std::string name, TFileDirectory *dir, Alignable *ali)=0
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)
virtual double	sumofweights ()=0
virtual int	type () const =0
int	useRes (int pattern=-1)
double	value (int parNum)
void	write (FILE *file, int which=0)
virtual	~MuonResidualsFitter ()

Protected Member Functions
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)
virtual void	inform (TMinuit *tMinuit)=0
void	initialize_table ()

Protected Attributes
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

2011/04/15 21:51:13

Revision

1.16

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

Id

MuonResidualsFitter.h,v 1.16 2011/04/15 21:51:13 khotilov Exp

Definition at line 84 of file MuonResidualsFitter.h.

Member Enumeration Documentation

anonymous enum

anonymous enum

Enumerator
kPureGaussian
kPowerLawTails
kROOTVoigt
kGaussPowerTails
kPureGaussian2D

Definition at line 86 of file MuonResidualsFitter.h.

{ kPureGaussian, kPowerLawTails, kROOTVoigt, kGaussPowerTails, kPureGaussian2D };

anonymous enum

anonymous enum

Enumerator
k1DOF
k5DOF
k6DOF
k6DOFrphi
kPositionFitter
kAngleFitter
kAngleBfieldFitter

Definition at line 88 of file MuonResidualsFitter.h.

{ k1DOF, k5DOF, k6DOF, k6DOFrphi, kPositionFitter, kAngleFitter, kAngleBfieldFitter };

anonymous enum

anonymous enum

Enumerator
k1111
k1110
k1100
k1010
k0010
k1000
k0100

Definition at line 90 of file MuonResidualsFitter.h.

{ k1111, k1110, k1100, k1010, k0010, k1000, k0100 };

Constructor & Destructor Documentation

MuonResidualsFitter()

MuonResidualsFitter::MuonResidualsFitter	(	int	residualsModel,
		int	minHits,
		int	useResiduals,
		bool	weightAlignment = true
	)

Definition at line 151 of file MuonResidualsFitter.cc.

References Exception, kGaussPowerTails, kPowerLawTails, kPureGaussian, kPureGaussian2D, kROOTVoigt, and m_residualsModel.

    : m_residualsModel(residualsModel),
      m_minHits(minHits),
      m_useResiduals(useResiduals),
      m_weightAlignment(weightAlignment),
      m_printLevel(0),
      m_strategy(1),
      m_cov(1),
      m_loglikelihood(0.) {
  if (m_residualsModel != kPureGaussian && m_residualsModel != kPowerLawTails && m_residualsModel != kROOTVoigt &&
      m_residualsModel != kGaussPowerTails && m_residualsModel != kPureGaussian2D)
    throw cms::Exception("MuonResidualsFitter") << "unrecognized residualsModel";
}

~MuonResidualsFitter()

MuonResidualsFitter::~MuonResidualsFitter ( )

virtual

Definition at line 165 of file MuonResidualsFitter.cc.

References residuals_begin(), and residuals_end().

                                          {
  for (std::vector<double *>::const_iterator residual = residuals_begin(); residual != residuals_end(); ++residual) {
    delete[](*residual);
  }
}

Member Function Documentation

computeHistogramRangeAndBinning()

void MuonResidualsFitter::computeHistogramRangeAndBinning	(	int	which,
		int &	nbins,
		double &	a,
		double &	b
	)

Definition at line 510 of file MuonResidualsFitter.cc.

References a, b, gather_cfg::cout, data, mps_fire::i, createfilelist::int, m_residuals, dqmiodumpmetadata::n, LaserClient_cfi::nbins, numResiduals(), funct::pow(), alignCSCRings::r, jetUpdater_cfi::sort, and eostools::which().

Referenced by histogramChi2GaussianFit().

                                                                                                     {
  // first, make a numeric array while discarding some crazy outliers
  double *data = new double[numResiduals()];
  int n = 0;
  for (std::vector<double *>::const_iterator r = m_residuals.begin(); r != m_residuals.end(); r++)
    if (fabs((*r)[which]) < 50.) {
      data[n] = (*r)[which];
      n++;
    }

  // compute "3 normal sigma" and regular interquantile ranges
  const int n_quantiles = 7;
  double probabilities[n_quantiles] = {0.00135, 0.02275, 0.25, 0.5, 0.75, 0.97725, 0.99865};  // "3 normal sigma"
  //double probabilities[n_quantiles] = {0.02275, 0.25, 0.75, 0.97725}; // "2 normal sigma"
  double quantiles[n_quantiles];
  std::sort(data, data + n);
  TMath::Quantiles(n, n_quantiles, data, quantiles, probabilities, true, nullptr, 7);
  delete[] data;
  double iqr = quantiles[4] - quantiles[2];

  // estimate optimal bin size according to Freedman-Diaconis rule
  double hbin = 2 * iqr / pow(n, 1. / 3);

  a = quantiles[1];
  b = quantiles[5];
  nbins = (int)((b - a) / hbin + 3.);  // add extra safety margin of 3

  std::cout << "   quantiles: ";
  for (int i = 0; i < n_quantiles; i++)
    std::cout << quantiles[i] << " ";
  std::cout << std::endl;
  //cout<<"n="<<select_count<<" quantiles ["<<quantiles[1]<<", "<<quantiles[2]<<"]  IQR="<<iqr
  //  <<"  full range=["<<minx<<","<<maxx<<"]"<<"  2 normal sigma quantile range = ["<<quantiles[0]<<", "<<quantiles[3]<<"]"<<endl;
  std::cout << "   optimal h=" << hbin << " nbins=" << nbins << std::endl;
}

correctBField() [1/2]

virtual void MuonResidualsFitter::correctBField ( )

pure virtual

Implemented in MuonResiduals6DOFFitter, MuonResiduals5DOFFitter, and MuonResiduals6DOFrphiFitter.

Referenced by MuonResiduals6DOFrphiFitter::correctBField(), MuonResiduals5DOFFitter::correctBField(), MuonResiduals6DOFFitter::correctBField(), and MuonResidualsTwoBin::correctBField().

correctBField() [2/2]

void MuonResidualsFitter::correctBField ( int  idx_momentum, int  idx_q  )

virtual

Definition at line 822 of file MuonResidualsFitter.cc.

References newFWLiteAna::bin, submitPVResolutionJobs::count, gather_cfg::cout, MillePedeFileConverter_cfg::e, mps_fire::i, createfilelist::int, dqmiolumiharvest::j, m_residuals, m_residuals_ok, cosmictrackSelector_cfi::min_pt, DiDispStaMuonMonitor_cfi::pt, alignCSCRings::r, residuals_begin(), residuals_end(), findQualityFiles::size, and jetUpdater_cfi::sort.

                                                                   {
  const int Nbin = 17;
  // find max 1/pt and bin width
  double min_pt = 9999.;
  for (std::vector<double *>::const_iterator r = residuals_begin(); r != residuals_end(); ++r) {
    double pt = fabs((*r)[idx_momentum]);
    if (pt < min_pt)
      min_pt = pt;
  }
  min_pt -= 0.01;  // to prevent bin # overflow
  const double bin_width = 1. / min_pt / Nbin;

  // fill indices of positive and negative charge residuals in each bin
  std::vector<size_t> pos[Nbin], neg[Nbin], to_erase;
  for (size_t i = 0; i < m_residuals.size(); i++) {
    if (!m_residuals_ok[i])
      continue;
    int bin = (int)floor(1. / fabs(m_residuals[i][idx_momentum]) / bin_width);
    if (m_residuals[i][idx_q] > 0)
      pos[bin].push_back(i);
    else
      neg[bin].push_back(i);
  }

  // equalize pos and neg in each bin
  for (int j = 0; j < Nbin; j++) {
    size_t psize = pos[j].size();
    size_t nsize = neg[j].size();
    if (psize == nsize)
      continue;

    std::set<int> idx_set;  // use a set to collect certain number of unique random indices to erase
    if (psize > nsize) {
      while (idx_set.size() < psize - nsize)
        idx_set.insert(gRandom->Integer(psize));
      for (std::set<int>::iterator it = idx_set.begin(); it != idx_set.end(); it++)
        to_erase.push_back(pos[j][*it]);
    } else {
      while (idx_set.size() < nsize - psize)
        idx_set.insert(gRandom->Integer(nsize));
      for (std::set<int>::iterator it = idx_set.begin(); it != idx_set.end(); it++)
        to_erase.push_back(neg[j][*it]);
    }
  }
  // sort in descending order, so we safely go from higher to lower indices:
  std::sort(to_erase.begin(), to_erase.end(), std::greater<int>());
  for (std::vector<size_t>::const_iterator e = to_erase.begin(); e != to_erase.end(); ++e) {
    m_residuals_ok[*e] = false;
    //delete[] *(m_residuals.begin() + *e);
    //m_residuals.erase(m_residuals.begin() + *e);
  }

  std::vector<size_t> apos[Nbin], aneg[Nbin];
  for (size_t i = 0; i < m_residuals.size(); i++) {
    if (!m_residuals_ok[i])
      continue;
    int bin = (int)floor(1. / fabs(m_residuals[i][idx_momentum]) / bin_width);
    if (m_residuals[i][idx_q] > 0)
      apos[bin].push_back(i);
    else
      aneg[bin].push_back(i);
  }
  for (int j = 0; j < Nbin; j++)
    std::cout << "bin " << j << ": [pos,neg] sizes before & after: [" << pos[j].size() << "," << neg[j].size()
              << "] -> [" << apos[j].size() << "," << aneg[j].size() << "]" << std::endl;
  std::cout << " N residuals " << m_residuals.size() << " -> "
            << (size_t)std::count(m_residuals_ok.begin(), m_residuals_ok.end(), true) << std::endl;
}

correlationMatrix()

TMatrixDSym MuonResidualsFitter::correlationMatrix

(

)

covarianceElement()

double MuonResidualsFitter::covarianceElement	(	int	parNum1,
		int	parNum2
	)

Definition at line 191 of file MuonResidualsFitter.cc.

References cms::cuda::assert(), m_cov, npar(), and parNum2parIdx().

                                                                      {
  assert(0 <= parNum1 && parNum1 < npar());
  assert(0 <= parNum2 && parNum2 < npar());
  assert(m_cov.GetNcols() == npar());  // m_cov might have not yet been resized to account for proper #parameters
  return m_cov(parNum2parIdx(parNum1), parNum2parIdx(parNum2));
}

covarianceMatrix()

TMatrixDSym MuonResidualsFitter::covarianceMatrix ( )

inline

Definition at line 158 of file MuonResidualsFitter.h.

References m_cov.

{ return m_cov; }

dofit()

bool MuonResidualsFitter::dofit ( void(*)(int &, double *, double &, double *, int)  fcn, 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  )

protected

Definition at line 272 of file MuonResidualsFitter.cc.

References MillePedeFileConverter_cfg::e, fcn(), fixed(), LaserClient_cfi::high, mps_fire::i, inform(), LaserClient_cfi::low, m_cov, m_error, m_loglikelihood, m_strategy, m_value, MuonResidualsFitter_TMinuit, npar(), and findQualityFiles::v.

Referenced by MuonResidualsAngleFitter::fit(), MuonResidualsPositionFitter::fit(), MuonResidualsBfieldAngleFitter::fit(), MuonResiduals1DOFFitter::fit(), MuonResiduals6DOFrphiFitter::fit(), MuonResiduals5DOFFitter::fit(), and MuonResiduals6DOFFitter::fit().

                                                         {
  MuonResidualsFitterFitInfo *fitinfo = new MuonResidualsFitterFitInfo(this);

  MuonResidualsFitter_TMinuit = new TMinuit(npar());
  // MuonResidualsFitter_TMinuit->SetPrintLevel(m_printLevel);
  MuonResidualsFitter_TMinuit->SetPrintLevel();
  MuonResidualsFitter_TMinuit->SetObjectFit(fitinfo);
  MuonResidualsFitter_TMinuit->SetFCN(fcn);
  inform(MuonResidualsFitter_TMinuit);

  std::vector<int>::const_iterator iNum = parNum.begin();
  std::vector<std::string>::const_iterator iName = parName.begin();
  std::vector<double>::const_iterator istart = start.begin();
  std::vector<double>::const_iterator istep = step.begin();
  std::vector<double>::const_iterator ilow = low.begin();
  std::vector<double>::const_iterator ihigh = high.begin();

  //MuonResidualsFitter_TMinuit->SetPrintLevel(-1);

  for (; iNum != parNum.end(); ++iNum, ++iName, ++istart, ++istep, ++ilow, ++ihigh) {
    MuonResidualsFitter_TMinuit->DefineParameter(*iNum, iName->c_str(), *istart, *istep, *ilow, *ihigh);
    if (fixed(*iNum))
      MuonResidualsFitter_TMinuit->FixParameter(*iNum);
  }

  double arglist[10];
  int ierflg;
  int smierflg;  //second MIGRAD ierflg

  // chi^2 errors should be 1.0, log-likelihood should be 0.5
  for (int i = 0; i < 10; i++)
    arglist[i] = 0.;
  arglist[0] = 0.5;
  ierflg = 0;
  smierflg = 0;
  MuonResidualsFitter_TMinuit->mnexcm("SET ERR", arglist, 1, ierflg);
  if (ierflg != 0) {
    delete MuonResidualsFitter_TMinuit;
    delete fitinfo;
    return false;
  }

  // set strategy = 2 (more refined fits)
  for (int i = 0; i < 10; i++)
    arglist[i] = 0.;
  arglist[0] = m_strategy;
  ierflg = 0;
  MuonResidualsFitter_TMinuit->mnexcm("SET STR", arglist, 1, ierflg);
  if (ierflg != 0) {
    delete MuonResidualsFitter_TMinuit;
    delete fitinfo;
    return false;
  }

  bool try_again = false;

  // minimize
  for (int i = 0; i < 10; i++)
    arglist[i] = 0.;
  arglist[0] = 50000;
  ierflg = 0;
  MuonResidualsFitter_TMinuit->mnexcm("MIGRAD", arglist, 1, ierflg);
  if (ierflg != 0)
    try_again = true;

  // just once more, if needed (using the final Minuit parameters from the failed fit; often works)
  if (try_again) {
    for (int i = 0; i < 10; i++)
      arglist[i] = 0.;
    arglist[0] = 50000;
    MuonResidualsFitter_TMinuit->mnexcm("MIGRAD", arglist, 1, smierflg);
  }

  Double_t fmin, fedm, errdef;
  Int_t npari, nparx, istat;
  MuonResidualsFitter_TMinuit->mnstat(fmin, fedm, errdef, npari, nparx, istat);

  if (istat != 3) {
    for (int i = 0; i < 10; i++)
      arglist[i] = 0.;
    ierflg = 0;
    MuonResidualsFitter_TMinuit->mnexcm("HESSE", arglist, 0, ierflg);
  }

  // read-out the results
  m_loglikelihood = -fmin;

  m_value.clear();
  m_error.clear();
  for (int i = 0; i < npar(); i++) {
    double v, e;
    MuonResidualsFitter_TMinuit->GetParameter(i, v, e);
    m_value.push_back(v);
    m_error.push_back(e);
  }
  m_cov.ResizeTo(npar(), npar());
  MuonResidualsFitter_TMinuit->mnemat(m_cov.GetMatrixArray(), npar());

  delete MuonResidualsFitter_TMinuit;
  delete fitinfo;
  if (smierflg != 0)
    return false;
  return true;
}

eraseNotSelectedResiduals()

void MuonResidualsFitter::eraseNotSelectedResiduals

(

)

Definition at line 891 of file MuonResidualsFitter.cc.

References submitPVResolutionJobs::count, gather_cfg::cout, mps_fire::i, m_residuals, m_residuals_ok, and createJobs::tmp.

Referenced by MuonResidualsTwoBin::eraseNotSelectedResiduals().

                                                    {
  // it should probably be faster then doing erase
  size_t n_ok = (size_t)std::count(m_residuals_ok.begin(), m_residuals_ok.end(), true);
  std::vector<double *> tmp(n_ok, nullptr);
  std::cout << "residuals sizes: all=" << m_residuals.size() << " good=" << n_ok << std::endl;
  int iok = 0;
  for (size_t i = 0; i < m_residuals.size(); i++) {
    if (!m_residuals_ok[i])
      continue;
    tmp[iok++] = m_residuals[i];
  }
  m_residuals.swap(tmp);

  std::vector<bool> tmp_ok(n_ok, true);
  m_residuals_ok.swap(tmp_ok);

  std::cout << "residuals size after eraseNotSelectedResiduals =" << m_residuals.size()
            << "  ok size=" << m_residuals_ok.size() << std::endl;
}

errorerror()

double MuonResidualsFitter::errorerror ( int  parNum )

inline

Definition at line 149 of file MuonResidualsFitter.h.

References cms::cuda::assert(), m_error, and npar().

Referenced by MuonResidualsTwoBin::errorerror(), MuonResiduals6DOFrphiFitter::plot(), MuonResiduals5DOFFitter::plot(), and MuonResiduals6DOFFitter::plot().

                                {
    assert(0 <= parNum && parNum < npar());
    return m_error[parNum];
  }

fiducialCuts()

void MuonResidualsFitter::fiducialCuts	(	double	xMin = -80.0,
		double	xMax = 80.0,
		double	yMin = -80.0,
		double	yMax = 80.0,
		bool	fidcut1 = false
	)

Definition at line 719 of file MuonResidualsFitter.cc.

References m_residuals_ok, alignCSCRings::r, residuals_begin(), residuals_end(), multiplicitycorr_cfi::xMax, photonAnalyzer_cfi::xMin, multiplicitycorr_cfi::yMax, and photonAnalyzer_cfi::yMin.

Referenced by MuonResidualsTwoBin::fiducialCuts().

                                                                                                       {
  int iResidual = -1;

  int n_station = 9999;
  int n_wheel = 9999;
  int n_sector = 9999;

  double positionX = 9999.;
  double positionY = 9999.;

  double chambw = 9999.;
  double chambl = 9999.;

  for (std::vector<double *>::const_iterator r = residuals_begin(); r != residuals_end(); ++r) {
    iResidual++;
    if (!m_residuals_ok[iResidual])
      continue;

    if ((*r)[15] >
        0.0001) {  // this value is greater than zero (chamber width) for 6DOFs stations 1,2,3 better to change for type()!!!
      n_station = (*r)[12];
      n_wheel = (*r)[13];
      n_sector = (*r)[14];
      positionX = (*r)[4];
      positionY = (*r)[5];
      chambw = (*r)[15];
      chambl = (*r)[16];
    } else {  // in case of 5DOF residual the residual object index is different
      n_station = (*r)[10];
      n_wheel = (*r)[11];
      n_sector = (*r)[12];
      positionX = (*r)[2];
      positionY = (*r)[3];
      chambw = (*r)[13];
      chambl = (*r)[14];
    }

    if (fidcut1) {  // this is the standard fiducial cut used so far 80x80 cm in x,y
      if (positionX >= xMax || positionX <= xMin)
        m_residuals_ok[iResidual] = false;
      if (positionY >= yMax || positionY <= yMin)
        m_residuals_ok[iResidual] = false;
    }

    // Implementation of new fiducial cut

    double dtrkchamx = (chambw / 2.) - positionX;  // variables to cut tracks on the edge of the chambers
    double dtrkchamy = (chambl / 2.) - positionY;

    if (!fidcut1) {
      if (n_station == 4) {
        if ((n_wheel == -1 && n_sector == 3) ||
            (n_wheel == 1 &&
             n_sector == 4)) {  // FOR SHORT CHAMBER LENGTH IN:  WHEEL 1 SECTOR 4  AND  WHEEL -1 SECTOR 3
          if ((n_sector == 1 || n_sector == 2 || n_sector == 3 || n_sector == 5 || n_sector == 6 || n_sector == 7 ||
               n_sector == 8 || n_sector == 12) &&
              ((dtrkchamx < 40 || dtrkchamx > 380) || (dtrkchamy < 40.0 || dtrkchamy > 170.0)))
            m_residuals_ok[iResidual] = false;
          if ((n_sector == 4 || n_sector == 13) &&
              ((dtrkchamx < 40 || dtrkchamx > 280) || (dtrkchamy < 40.0 || dtrkchamy > 170.0)))
            m_residuals_ok[iResidual] = false;
          if ((n_sector == 9 || n_sector == 11) &&
              ((dtrkchamx < 40 || dtrkchamx > 180) || (dtrkchamy < 40.0 || dtrkchamy > 170.0)))
            m_residuals_ok[iResidual] = false;
          if ((n_sector == 10 || n_sector == 14) &&
              ((dtrkchamx < 40 || dtrkchamx > 220) || (dtrkchamy < 40.0 || dtrkchamy > 170.0)))
            m_residuals_ok[iResidual] = false;
        } else {
          if ((n_sector == 1 || n_sector == 2 || n_sector == 3 || n_sector == 5 || n_sector == 6 || n_sector == 7 ||
               n_sector == 8 || n_sector == 12) &&
              ((dtrkchamx < 40 || dtrkchamx > 380) || (dtrkchamy < 40.0 || dtrkchamy > 210.0)))
            m_residuals_ok[iResidual] = false;
          if ((n_sector == 4 || n_sector == 13) &&
              ((dtrkchamx < 40 || dtrkchamx > 280) || (dtrkchamy < 40.0 || dtrkchamy > 210.0)))
            m_residuals_ok[iResidual] = false;
          if ((n_sector == 9 || n_sector == 11) &&
              ((dtrkchamx < 40 || dtrkchamx > 180) || (dtrkchamy < 40.0 || dtrkchamy > 210.0)))
            m_residuals_ok[iResidual] = false;
          if ((n_sector == 10 || n_sector == 14) &&
              ((dtrkchamx < 40 || dtrkchamx > 220) || (dtrkchamy < 40.0 || dtrkchamy > 210.0)))
            m_residuals_ok[iResidual] = false;
        }
      } else {
        if ((n_wheel == -1 && n_sector == 3) || (n_wheel == 1 && n_sector == 4)) {
          if (n_station == 1 && ((dtrkchamx < 30.0 || dtrkchamx > 190.0) || (dtrkchamy < 40.0 || dtrkchamy > 170.0)))
            m_residuals_ok[iResidual] = false;
          if (n_station == 2 && ((dtrkchamx < 30.0 || dtrkchamx > 240.0) || (dtrkchamy < 40.0 || dtrkchamy > 170.0)))
            m_residuals_ok[iResidual] = false;
          if (n_station == 3 && ((dtrkchamx < 30.0 || dtrkchamx > 280.0) || (dtrkchamy < 40.0 || dtrkchamy > 170.0)))
            m_residuals_ok[iResidual] = false;
        } else {
          if (n_station == 1 && ((dtrkchamx < 30.0 || dtrkchamx > 190.0) || (dtrkchamy < 40.0 || dtrkchamy > 210.0)))
            m_residuals_ok[iResidual] = false;
          if (n_station == 2 && ((dtrkchamx < 30.0 || dtrkchamx > 240.0) || (dtrkchamy < 40.0 || dtrkchamy > 210.0)))
            m_residuals_ok[iResidual] = false;
          if (n_station == 3 && ((dtrkchamx < 30.0 || dtrkchamx > 280.0) || (dtrkchamy < 40.0 || dtrkchamy > 210.0)))
            m_residuals_ok[iResidual] = false;
        }
      }
    }
  }
}

fill()

void MuonResidualsFitter::fill

(

double *

residual

)

Definition at line 186 of file MuonResidualsFitter.cc.

References m_residuals, and m_residuals_ok.

Referenced by MuonResidualsTwoBin::fill(), read(), MuonResiduals6DOFrphiFitter::readNtuple(), MuonResiduals5DOFFitter::readNtuple(), and MuonResiduals6DOFFitter::readNtuple().

                                               {
  m_residuals.push_back(residual);
  m_residuals_ok.push_back(true);
}

fit()

virtual bool MuonResidualsFitter::fit ( Alignable *  ali )

pure virtual

Implemented in MuonResiduals6DOFFitter, MuonResiduals5DOFFitter, MuonResiduals6DOFrphiFitter, MuonResiduals1DOFFitter, MuonResidualsAngleFitter, MuonResidualsBfieldAngleFitter, and MuonResidualsPositionFitter.

Referenced by MuonResidualsTwoBin::fit(), and trackingPlots.Iteration::modules().

fix()

void MuonResidualsFitter::fix	(	int	parNum,
		bool	dofix = true
	)

Definition at line 171 of file MuonResidualsFitter.cc.

References cms::cuda::assert(), m_fixed, and npar().

Referenced by MuonResiduals6DOFrphiFitter::fit(), MuonResiduals5DOFFitter::fit(), MuonResiduals6DOFFitter::fit(), and MuonResidualsTwoBin::fix().

                                                    {
  assert(0 <= parNum && parNum < npar());
  if (m_fixed.empty())
    m_fixed.resize(npar(), false);
  m_fixed[parNum] = dofix;
}

fixed()

bool MuonResidualsFitter::fixed

(

int

parNum

)

Definition at line 178 of file MuonResidualsFitter.cc.

References cms::cuda::assert(), m_fixed, and npar().

Referenced by dofit(), MuonResiduals1DOFFitter::fit(), and MuonResidualsTwoBin::fixed().

                                          {
  assert(0 <= parNum && parNum < npar());
  if (m_fixed.empty())
    return false;
  else
    return m_fixed[parNum];
}

histogramChi2GaussianFit()

void MuonResidualsFitter::histogramChi2GaussianFit	(	int	which,
		double &	fit_mean,
		double &	fit_sigma
	)

Definition at line 546 of file MuonResidualsFitter.cc.

References a, b, computeHistogramRangeAndBinning(), gather_cfg::cout, DeadROC_duringRun::f1, gpuVertexFinder::hist, m_residuals, LaserClient_cfi::nbins, alignCSCRings::r, and eostools::which().

Referenced by selectPeakResiduals().

                                                                                                 {
  int nbins;
  double a, b;
  computeHistogramRangeAndBinning(which, nbins, a, b);
  if (a == b || a > b) {
    fit_mean = a;
    fit_sigma = 0;
    return;
  }

  TH1D *hist = new TH1D("htmp", "", nbins, a, b);
  for (std::vector<double *>::const_iterator r = m_residuals.begin(); r != m_residuals.end(); ++r)
    hist->Fill((*r)[which]);

  // do simple chi2 gaussian fit
  TF1 *f1 = new TF1("f1", "gaus", a, b);
  f1->SetParameter(0, hist->GetEntries());
  f1->SetParameter(1, 0);
  f1->SetParameter(2, hist->GetRMS());
  hist->Fit("f1", "RQ");
  // hist->Fit(f1,"RQ");

  fit_mean = f1->GetParameter(1);
  fit_sigma = f1->GetParameter(2);
  std::cout << " h(" << nbins << "," << a << "," << b << ") mu=" << fit_mean << " sig=" << fit_sigma << std::endl;

  delete f1;
  delete hist;
}

inform()

virtual void MuonResidualsFitter::inform ( TMinuit *  tMinuit )

protectedpure virtual

Implemented in MuonResiduals6DOFFitter, MuonResiduals5DOFFitter, MuonResiduals6DOFrphiFitter, MuonResiduals1DOFFitter, MuonResidualsAngleFitter, MuonResidualsBfieldAngleFitter, and MuonResidualsPositionFitter.

Referenced by dofit().

initialize_table()

void MuonResidualsFitter::initialize_table ( )

protected

Definition at line 205 of file MuonResidualsFitter.cc.

References gather_cfg::cout, Exception, kPowerLawTails, MuonResidualsFitter_compute_log_convolution(), MuonResidualsFitter_gsbinsize, MuonResidualsFitter_lookup_table, MuonResidualsFitter_numgsbins, MuonResidualsFitter_numtsbins, MuonResidualsFitter_table_initialized, MuonResidualsFitter_tsbinsize, residualsModel(), and heppy_batch::val.

Referenced by MuonResidualsAngleFitter::fit(), MuonResidualsPositionFitter::fit(), MuonResidualsBfieldAngleFitter::fit(), MuonResiduals1DOFFitter::fit(), MuonResiduals6DOFrphiFitter::fit(), MuonResiduals5DOFFitter::fit(), and MuonResiduals6DOFFitter::fit().

                                           {
  if (MuonResidualsFitter_table_initialized || residualsModel() != kPowerLawTails)
    return;
  MuonResidualsFitter_table_initialized = true;

  std::ifstream convolution_table("convolution_table.txt");
  if (convolution_table.is_open()) {
    int numgsbins = 0;
    int numtsbins = 0;
    double tsbinsize = 0.;
    double gsbinsize = 0.;

    convolution_table >> numgsbins >> numtsbins >> tsbinsize >> gsbinsize;
    if (numgsbins != MuonResidualsFitter_numgsbins || numtsbins != MuonResidualsFitter_numtsbins ||
        tsbinsize != MuonResidualsFitter_tsbinsize || gsbinsize != MuonResidualsFitter_gsbinsize) {
      throw cms::Exception("MuonResidualsFitter") << "convolution_table.txt has the wrong bin width/bin size.  Throw "
                                                     "it away and let the fitter re-create the file.\n";
    }

    for (int gsbin = 0; gsbin < MuonResidualsFitter_numgsbins; gsbin++) {
      for (int tsbin = 0; tsbin < MuonResidualsFitter_numtsbins; tsbin++) {
        int read_gsbin = 0;
        int read_tsbin = 0;
        double val = 0.;

        convolution_table >> read_gsbin >> read_tsbin >> val;
        if (read_gsbin != gsbin || read_tsbin != tsbin) {
          throw cms::Exception("MuonResidualsFitter")
              << "convolution_table.txt is out of order.  Throw it away and let the fitter re-create the file.\n";
        }
        MuonResidualsFitter_lookup_table[gsbin][tsbin] = val;
      }
    }
    convolution_table.close();
  } else {
    std::ofstream convolution_table2("convolution_table.txt");

    if (!convolution_table2.is_open())
      throw cms::Exception("MuonResidualsFitter") << "Couldn't write to file convolution_table.txt\n";

    convolution_table2 << MuonResidualsFitter_numgsbins << " " << MuonResidualsFitter_numtsbins << " "
                       << MuonResidualsFitter_tsbinsize << " " << MuonResidualsFitter_gsbinsize << std::endl;

    std::cout << "Initializing convolution look-up table (takes a few minutes)..." << std::endl;

    for (int gsbin = 0; gsbin < MuonResidualsFitter_numgsbins; gsbin++) {
      double gammaoversigma = double(gsbin) * MuonResidualsFitter_gsbinsize;
      std::cout << "    gsbin " << gsbin << "/" << MuonResidualsFitter_numgsbins << std::endl;
      for (int tsbin = 0; tsbin < MuonResidualsFitter_numtsbins; tsbin++) {
        double toversigma = double(tsbin) * MuonResidualsFitter_tsbinsize;

        // 1e-6 errors (out of a value of ~0.01) with max=100, step=0.001, power=4 (max=1000 does a little better with the tails)
        MuonResidualsFitter_lookup_table[gsbin][tsbin] =
            MuonResidualsFitter_compute_log_convolution(toversigma, gammaoversigma);

        // <10% errors with max=20, step=0.005, power=4 (faster computation for testing)
        // MuonResidualsFitter_lookup_table[gsbin][tsbin] = MuonResidualsFitter_compute_log_convolution(toversigma, gammaoversigma, 100., 0.005, 4.);

        convolution_table2 << gsbin << " " << tsbin << " " << MuonResidualsFitter_lookup_table[gsbin][tsbin]
                           << std::endl;
      }
    }
    convolution_table2.close();
    std::cout << "Initialization done!" << std::endl;
  }
}

loglikelihood()

double MuonResidualsFitter::loglikelihood ( )

inline

Definition at line 162 of file MuonResidualsFitter.h.

References m_loglikelihood.

Referenced by MuonResidualsTwoBin::loglikelihood().

{ return m_loglikelihood; }

ndata()

virtual int MuonResidualsFitter::ndata ( )

pure virtual

Implemented in MuonResiduals6DOFFitter, MuonResiduals5DOFFitter, MuonResiduals6DOFrphiFitter, MuonResiduals1DOFFitter, MuonResidualsAngleFitter, MuonResidualsBfieldAngleFitter, and MuonResidualsPositionFitter.

Referenced by MuonResidualsTwoBin::ndata(), read(), and write().

nfixed()

int MuonResidualsFitter::nfixed ( )

inline

Definition at line 130 of file MuonResidualsFitter.h.

References submitPVResolutionJobs::count, and m_fixed.

{ return std::count(m_fixed.begin(), m_fixed.end(), true); }

npar()

virtual int MuonResidualsFitter::npar ( )

pure virtual

Implemented in MuonResiduals6DOFFitter, MuonResiduals5DOFFitter, MuonResiduals6DOFrphiFitter, MuonResiduals1DOFFitter, MuonResidualsAngleFitter, MuonResidualsBfieldAngleFitter, and MuonResidualsPositionFitter.

Referenced by covarianceElement(), dofit(), errorerror(), fix(), fixed(), MuonResidualsTwoBin::npar(), and value().

numResiduals()

long MuonResidualsFitter::numResiduals ( ) const

inline

Definition at line 125 of file MuonResidualsFitter.h.

References m_residuals.

Referenced by computeHistogramRangeAndBinning(), MuonResidualsTwoBin::numResidualsNeg(), MuonResidualsTwoBin::numResidualsPos(), MuonResidualsAngleFitter::plot(), MuonResidualsPositionFitter::plot(), MuonResidualsBfieldAngleFitter::plot(), read(), selectPeakResiduals(), selectPeakResiduals_simple(), MuonResidualsAngleFitter::sumofweights(), MuonResidualsPositionFitter::sumofweights(), MuonResidualsBfieldAngleFitter::sumofweights(), and write().

{ return m_residuals.size(); }

numsegments()

long MuonResidualsFitter::numsegments

(

)

Definition at line 198 of file MuonResidualsFitter.cc.

References EgammaValidation_cff::num, residuals_begin(), and residuals_end().

Referenced by MuonResidualsTwoBin::numsegments().

                                      {
  long num = 0;
  for (std::vector<double *>::const_iterator resiter = residuals_begin(); resiter != residuals_end(); ++resiter)
    num++;
  return num;
}

parNum2parIdx()

int MuonResidualsFitter::parNum2parIdx ( int  parNum )

inline

Definition at line 156 of file MuonResidualsFitter.h.

References m_parNum2parIdx.

Referenced by covarianceElement().

{ return m_parNum2parIdx[parNum]; }

plot()

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

pure virtual

Implemented in MuonResiduals6DOFFitter, MuonResiduals5DOFFitter, MuonResiduals6DOFrphiFitter, MuonResiduals1DOFFitter, MuonResidualsAngleFitter, MuonResidualsBfieldAngleFitter, and MuonResidualsPositionFitter.

Referenced by MuonResidualsTwoBin::plot().

plotsimple()

void MuonResidualsFitter::plotsimple	(	std::string	name,
		TFileDirectory *	dir,
		int	which,
		double	multiplier
	)

Definition at line 474 of file MuonResidualsFitter.cc.

References DeadROC_duringRun::dir, gpuVertexFinder::hist, Skims_PA_cff::name, alignCSCRings::r, residuals_begin(), residuals_end(), eostools::which(), and svgfig::window().

Referenced by MuonResidualsTwoBin::plotsimple().

                                                                                                      {
  double window = 100.;
  if (which == 0)
    window = 2. * 30.;
  else if (which == 1)
    window = 2. * 30.;
  else if (which == 2)
    window = 2. * 20.;
  else if (which == 3)
    window = 2. * 50.;

  TH1F *hist = dir->make<TH1F>(name.c_str(), "", 200, -window, window);
  for (std::vector<double *>::const_iterator r = residuals_begin(); r != residuals_end(); ++r)
    hist->Fill(multiplier * (*r)[which]);
}

plotweighted()

void MuonResidualsFitter::plotweighted	(	std::string	name,
		TFileDirectory *	dir,
		int	which,
		int	whichredchi2,
		double	multiplier
	)

Definition at line 490 of file MuonResidualsFitter.cc.

References DeadROC_duringRun::dir, gpuVertexFinder::hist, Skims_PA_cff::name, alignCSCRings::r, residuals_begin(), residuals_end(), eostools::which(), and svgfig::window().

Referenced by MuonResidualsTwoBin::plotweighted().

                                                                                         {
  double window = 100.;
  if (which == 0)
    window = 2. * 30.;
  else if (which == 1)
    window = 2. * 30.;
  else if (which == 2)
    window = 2. * 20.;
  else if (which == 3)
    window = 2. * 50.;

  TH1F *hist = dir->make<TH1F>(name.c_str(), "", 200, -window, window);
  for (std::vector<double *>::const_iterator r = residuals_begin(); r != residuals_end(); ++r) {
    double weight = 1. / (*r)[whichredchi2];
    if (TMath::Prob(1. / weight * 12, 12) < 0.99)
      hist->Fill(multiplier * (*r)[which], weight);
  }
}

read()

void MuonResidualsFitter::read	(	FILE *	file,
		int	which = 0
	)

Definition at line 418 of file MuonResidualsFitter.cc.

References MillePedeFileConverter_cfg::e, Exception, geometryDiff::file, fill(), mps_fire::i, m_residuals, m_residuals_ok, ndata(), numResiduals(), postprocess-scan-build::rows, and eostools::which().

Referenced by edmIntegrityCheck.PublishToFileSystem::get(), and MuonResidualsTwoBin::read().

                                                    {
  long rows = -100;
  int cols = -100;
  int readwhich = -100;

  fread(&rows, sizeof(long), 1, file);
  fread(&cols, sizeof(int), 1, file);
  fread(&readwhich, sizeof(int), 1, file);

  if (cols != ndata() || rows < 0 || readwhich != which) {
    throw cms::Exception("MuonResidualsFitter")
        << "temporary file is corrupted (which = " << which << " readwhich = " << readwhich << " rows = " << rows
        << " cols = " << cols << ")\n";
  }

  double *likeAChecksum = new double[cols];
  double *likeAChecksum2 = new double[cols];
  for (int i = 0; i < cols; i++) {
    likeAChecksum[i] = 0.;
    likeAChecksum2[i] = 0.;
  }

  m_residuals.reserve(rows);
  for (long row = 0; row < rows; row++) {
    double *residual = new double[cols];
    fread(residual, sizeof(double), cols, file);
    fill(residual);
    for (int i = 0; i < cols; i++) {
      if (fabs(residual[i]) > likeAChecksum[i])
        likeAChecksum[i] = fabs(residual[i]);
      if (fabs(residual[i]) < likeAChecksum2[i])
        likeAChecksum2[i] = fabs(residual[i]);
    }
  }  // end loop over records in file

  double *readChecksum = new double[cols];
  double *readChecksum2 = new double[cols];
  fread(readChecksum, sizeof(double), cols, file);
  fread(readChecksum2, sizeof(double), cols, file);

  for (int i = 0; i < cols; i++) {
    if (fabs(likeAChecksum[i] - readChecksum[i]) > 1e-10 ||
        fabs(1. / likeAChecksum2[i] - 1. / readChecksum2[i]) > 1e10) {
      throw cms::Exception("MuonResidualsFitter")
          << "temporary file is corrupted (which = " << which << " rows = " << rows << " likeAChecksum "
          << likeAChecksum[i] << " != readChecksum " << readChecksum[i] << " "
          << " likeAChecksum2 " << likeAChecksum2[i] << " != readChecksum2 " << readChecksum2[i] << ")\n";
    }
  }

  m_residuals_ok.resize(numResiduals(), true);

  delete[] likeAChecksum;
  delete[] likeAChecksum2;
}

residuals_begin()

std::vector<double *>::const_iterator MuonResidualsFitter::residuals_begin ( ) const

inline

Definition at line 190 of file MuonResidualsFitter.h.

References m_residuals.

Referenced by correctBField(), fiducialCuts(), MuonResidualsBfieldAngleFitter::fit(), MuonResidualsAngleFitter::fit(), MuonResidualsPositionFitter::fit(), MuonResiduals1DOFFitter::fit(), MuonResiduals1DOFFitter_FCN(), MuonResiduals5DOFFitter_FCN(), MuonResiduals6DOFFitter_FCN(), MuonResiduals6DOFrphiFitter_FCN(), MuonResidualsAngleFitter_FCN(), MuonResidualsBfieldAngleFitter_FCN(), MuonResidualsPositionFitter_FCN(), numsegments(), MuonResidualsAngleFitter::plot(), MuonResidualsBfieldAngleFitter::plot(), MuonResiduals1DOFFitter::plot(), MuonResidualsPositionFitter::plot(), MuonResiduals6DOFrphiFitter::plot(), MuonResiduals5DOFFitter::plot(), MuonResiduals6DOFFitter::plot(), plotsimple(), plotweighted(), MuonResidualsTwoBin::residualsNeg_begin(), MuonResidualsTwoBin::residualsPos_begin(), MuonResiduals1DOFFitter::sumofweights(), MuonResiduals6DOFrphiFitter::sumofweights(), MuonResiduals5DOFFitter::sumofweights(), MuonResiduals6DOFFitter::sumofweights(), write(), and ~MuonResidualsFitter().

{ return m_residuals.begin(); }

residuals_end()

std::vector<double *>::const_iterator MuonResidualsFitter::residuals_end ( ) const

inline

Definition at line 191 of file MuonResidualsFitter.h.

References m_residuals.

Referenced by correctBField(), fiducialCuts(), MuonResidualsBfieldAngleFitter::fit(), MuonResidualsAngleFitter::fit(), MuonResidualsPositionFitter::fit(), MuonResiduals1DOFFitter::fit(), MuonResiduals1DOFFitter_FCN(), MuonResiduals5DOFFitter_FCN(), MuonResiduals6DOFFitter_FCN(), MuonResiduals6DOFrphiFitter_FCN(), MuonResidualsAngleFitter_FCN(), MuonResidualsBfieldAngleFitter_FCN(), MuonResidualsPositionFitter_FCN(), numsegments(), MuonResidualsAngleFitter::plot(), MuonResidualsBfieldAngleFitter::plot(), MuonResiduals1DOFFitter::plot(), MuonResidualsPositionFitter::plot(), MuonResiduals6DOFrphiFitter::plot(), MuonResiduals5DOFFitter::plot(), MuonResiduals6DOFFitter::plot(), plotsimple(), plotweighted(), MuonResidualsTwoBin::residualsNeg_end(), MuonResidualsTwoBin::residualsPos_end(), MuonResiduals1DOFFitter::sumofweights(), MuonResiduals6DOFrphiFitter::sumofweights(), MuonResiduals5DOFFitter::sumofweights(), MuonResiduals6DOFFitter::sumofweights(), write(), and ~MuonResidualsFitter().

{ return m_residuals.end(); }

residualsModel()

int MuonResidualsFitter::residualsModel ( ) const

inline

Definition at line 124 of file MuonResidualsFitter.h.

References m_residualsModel.

Referenced by MuonResidualsAngleFitter::fit(), MuonResidualsBfieldAngleFitter::fit(), MuonResidualsPositionFitter::fit(), MuonResiduals1DOFFitter::fit(), MuonResiduals6DOFrphiFitter::fit(), MuonResiduals5DOFFitter::fit(), MuonResiduals6DOFFitter::fit(), initialize_table(), MuonResiduals1DOFFitter_FCN(), MuonResiduals5DOFFitter_FCN(), MuonResiduals6DOFFitter_FCN(), MuonResiduals6DOFrphiFitter_FCN(), MuonResidualsAngleFitter_FCN(), MuonResidualsBfieldAngleFitter_FCN(), MuonResidualsPositionFitter_FCN(), MuonResidualsAngleFitter::npar(), MuonResiduals1DOFFitter::npar(), MuonResidualsBfieldAngleFitter::npar(), MuonResidualsPositionFitter::npar(), MuonResiduals6DOFrphiFitter::npar(), MuonResiduals5DOFFitter::npar(), MuonResiduals6DOFFitter::npar(), MuonResidualsPositionFitter::plot(), MuonResidualsBfieldAngleFitter::plot(), MuonResidualsAngleFitter::plot(), MuonResiduals1DOFFitter::plot(), MuonResiduals6DOFrphiFitter::plot(), MuonResiduals5DOFFitter::plot(), MuonResiduals6DOFFitter::plot(), and MuonResidualsTwoBin::residualsModel().

{ return m_residualsModel; }

selectedResidualsFlags()

std::vector<bool>& MuonResidualsFitter::selectedResidualsFlags ( )

inline

Definition at line 209 of file MuonResidualsFitter.h.

References m_residuals_ok.

Referenced by MuonResidualsTwoBin::residualsNeg_ok_begin(), MuonResidualsTwoBin::residualsNeg_ok_end(), MuonResidualsTwoBin::residualsPos_ok_begin(), and MuonResidualsTwoBin::residualsPos_ok_end().

{ return m_residuals_ok; }

selectPeakResiduals()

void MuonResidualsFitter::selectPeakResiduals	(	double	nsigma,
		int	nvar,
		int *	vars
	)

Definition at line 579 of file MuonResidualsFitter.cc.

References cms::cuda::assert(), gather_cfg::cout, histogramChi2GaussianFit(), m_center, m_radii, m_residuals, m_residuals_ok, numResiduals(), funct::pow(), alignCSCRings::r, findQualityFiles::v, and eostools::which().

Referenced by MuonResidualsTwoBin::selectPeakResiduals().

                                                                                {
  // does not make sense for small statistics
  if (numResiduals() < 25)
    return;

  int nbefore = numResiduals();

  //just to be sure (can't see why it might ever be more then 10)
  assert(nvar <= 10);

  // estimate nvar-D ellipsoid center & axes
  for (int v = 0; v < nvar; v++) {
    int which = vars[v];
    histogramChi2GaussianFit(which, m_center[which], m_radii[which]);
    m_radii[which] = nsigma * m_radii[which];
  }

  // filter out residuals that don't fit into the ellipsoid
  std::vector<double *>::iterator r = m_residuals.begin();
  while (r != m_residuals.end()) {
    double ellipsoid_sum = 0;
    for (int v = 0; v < nvar; v++) {
      int which = vars[v];
      if (m_radii[which] == 0.)
        continue;
      ellipsoid_sum += pow(((*r)[which] - m_center[which]) / m_radii[which], 2);
    }
    if (ellipsoid_sum <= 1.)
      ++r;
    else {
      m_residuals_ok[r - m_residuals.begin()] = false;
      ++r;
      // delete [] (*r);
      // r = m_residuals.erase(r);
    }
  }
  std::cout << " N residuals " << nbefore << " -> " << numResiduals() << std::endl;
}

selectPeakResiduals_simple()

void MuonResidualsFitter::selectPeakResiduals_simple	(	double	nsigma,
		int	nvar,
		int *	vars
	)

Definition at line 620 of file MuonResidualsFitter.cc.

References cms::cuda::assert(), submitPVResolutionJobs::count, gather_cfg::cout, data, HLT_2022v12_cff::distance, mps_fire::i, m_residuals, m_residuals_ok, numResiduals(), alignCSCRings::r, mathSSE::sqrt(), and findQualityFiles::v.

                                                                                       {
  //std::cout<<"doing selectpeakresiduals: nsig="<<nsigma<<" nvar="<<nvar<<" vars=";
  for (int i = 0; i < nvar; ++i)
    std::cout << vars[i] << " ";
  std::cout << std::endl;

  // does not make sense for small statistics set to 50
  // YP changed it to 10 for test
  if (numResiduals() < 10)
    return;
  // if (numResiduals()<50) return;

  size_t nbefore = numResiduals();
  std::cout << " N residuals " << nbefore << " ~ "
            << (size_t)std::count(m_residuals_ok.begin(), m_residuals_ok.end(), true) << std::endl;
  //just to be sure (can't see why it might ever be more then 10)
  assert(nvar <= 10 && nvar > 0);

  std::vector<double *>::iterator r = m_residuals.begin();

  // it's awkward, but the 1D case has to be handled separately
  if (nvar == 1) {
    std::cout << "1D case" << std::endl;
    // get robust estimates for the peak and sigma
    double *data = new double[nbefore];
    for (size_t i = 0; i < nbefore; i++)
      data[i] = m_residuals[i][vars[0]];
    double peak, sigma;
    TRobustEstimator re;
    re.EvaluateUni(nbefore, data, peak, sigma);

    // filter out residuals that are more then nsigma away from the peak
    while (r != m_residuals.end()) {
      double distance = fabs(((*r)[vars[0]] - peak) / sigma);
      if (distance <= nsigma)
        ++r;
      else {
        m_residuals_ok[r - m_residuals.begin()] = false;
        ++r;
        //delete [] (*r);
        //r = m_residuals.erase(r);
      }
    }
    std::cout << " N residuals " << nbefore << " -> " << numResiduals() << std::endl;
    return;
  }  // end 1D case

  // initialize and run the robust estimator for D>1
  std::cout << "D>1 case" << std::endl;
  TRobustEstimator re(nbefore + 1, nvar);
  std::cout << "nbefore " << nbefore << " nvar " << nvar << std::endl;
  r = m_residuals.begin();
  std::cout << "+++++ JUST before loop while (r != m_residuals.end())" << std::endl;
  int counter1 = 0;
  while (r != m_residuals.end()) {
    double *row = new double[nvar];
    for (int v = 0; v < nvar; v++)
      row[v] = (*r)[vars[v]];
    re.AddRow(row);
    // delete[] row;
    ++r;
    counter1++;
  }
  std::cout << "counter1 " << counter1 << std::endl;
  std::cout << "+++++ JUST after loop while (r != m_residuals.end())" << std::endl;
  re.Evaluate();
  std::cout << "+++++ JUST after re.Evaluate()" << std::endl;

  // get nvar-dimensional ellipsoid center & covariance
  TVectorD M(nvar);
  re.GetMean(M);
  TMatrixDSym Cov(nvar);
  re.GetCovariance(Cov);
  Cov.Invert();

  // calculate the normalized radius for this nvar-dimensional ellipsoid from a 1D-Gaussian nsigma equivalent distance
  double conf_1d = TMath::Erf(nsigma / sqrt(2));
  double surf_radius = sqrt(TMath::ChisquareQuantile(conf_1d, nvar));

  // filter out residuals that are outside of the covariance ellipsoid with the above normalized radius
  r = m_residuals.begin();
  while (r != m_residuals.end()) {
    TVectorD res(nvar);
    for (int v = 0; v < nvar; v++)
      res[v] = (*r)[vars[v]];
    double distance = sqrt(Cov.Similarity(res - M));
    if (distance <= surf_radius)
      ++r;
    else {
      m_residuals_ok[r - m_residuals.begin()] = false;
      ++r;
      //delete [] (*r);
      //r = m_residuals.erase(r);
    }
  }
  std::cout << " N residuals " << nbefore << " -> "
            << (size_t)std::count(m_residuals_ok.begin(), m_residuals_ok.end(), true) << std::endl;
}

setInitialValue()

void MuonResidualsFitter::setInitialValue ( int  parNum, double  value  )

inline

Definition at line 135 of file MuonResidualsFitter.h.

References m_parNum2InitValue, and value().

{ m_parNum2InitValue[parNum] = value; }

setPrintLevel()

void MuonResidualsFitter::setPrintLevel ( int  printLevel )

inline

Definition at line 132 of file MuonResidualsFitter.h.

References m_printLevel.

Referenced by MuonResidualsTwoBin::setPrintLevel().

{ m_printLevel = printLevel; }

setStrategy()

void MuonResidualsFitter::setStrategy ( int  strategy )

inline

Definition at line 133 of file MuonResidualsFitter.h.

References m_strategy.

Referenced by MuonResidualsTwoBin::setStrategy().

{ m_strategy = strategy; }

sumofweights()

virtual double MuonResidualsFitter::sumofweights ( )

pure virtual

Implemented in MuonResiduals6DOFFitter, MuonResiduals5DOFFitter, MuonResiduals6DOFrphiFitter, MuonResidualsAngleFitter, MuonResidualsBfieldAngleFitter, MuonResidualsPositionFitter, and MuonResiduals1DOFFitter.

Referenced by MuonResidualsTwoBin::sumofweights().

type()

virtual int MuonResidualsFitter::type ( ) const

pure virtual

Implemented in MuonResiduals6DOFFitter, MuonResiduals5DOFFitter, MuonResiduals6DOFrphiFitter, MuonResiduals1DOFFitter, MuonResidualsAngleFitter, MuonResidualsBfieldAngleFitter, and MuonResidualsPositionFitter.

Referenced by MuonResidualsTwoBin::type().

useRes()

int MuonResidualsFitter::useRes ( int  pattern = -1 )

inline

Definition at line 119 of file MuonResidualsFitter.h.

References m_useResiduals, and topSingleLeptonDQM_PU_cfi::pattern.

Referenced by MuonResiduals6DOFrphiFitter::fit(), MuonResiduals5DOFFitter::fit(), MuonResiduals6DOFFitter::fit(), MuonResiduals5DOFFitter_FCN(), MuonResiduals6DOFFitter_FCN(), MuonResiduals6DOFrphiFitter_FCN(), and MuonResidualsTwoBin::useRes().

                               {
    if (pattern >= 0)
      m_useResiduals = pattern;
    return m_useResiduals;
  }

value()

double MuonResidualsFitter::value ( int  parNum )

inline

Definition at line 145 of file MuonResidualsFitter.h.

References cms::cuda::assert(), m_value, and npar().

Referenced by Types.int32::__nonzero__(), Types.uint32::__nonzero__(), Types.int64::__nonzero__(), Types.uint64::__nonzero__(), Types.double::__nonzero__(), Types.bool::__nonzero__(), Types.string::__nonzero__(), MuonResidualsTwoBin::antisym(), average.Average::average(), Types.string::configValue(), Types.FileInPath::configValue(), Mixins.UsingBlock::dumpPython(), Mixins.UsingBlock::insertInto(), Types.int32::insertInto(), Types.uint32::insertInto(), Types.int64::insertInto(), Types.uint64::insertInto(), Types.double::insertInto(), Types.bool::insertInto(), Types.string::insertInto(), Types.FileInPath::insertInto(), Types.vint32::insertInto(), Types.vuint32::insertInto(), Types.vint64::insertInto(), Types.vuint64::insertInto(), Types.vdouble::insertInto(), Types.vbool::insertInto(), Types.vstring::insertInto(), MuonResidualsPositionFitter::plot(), MuonResidualsBfieldAngleFitter::plot(), MuonResiduals1DOFFitter::plot(), MuonResidualsAngleFitter::plot(), MuonResiduals6DOFrphiFitter::plot(), MuonResiduals5DOFFitter::plot(), MuonResiduals6DOFFitter::plot(), setInitialValue(), and MuonResidualsTwoBin::value().

                           {
    assert(0 <= parNum && parNum < npar());
    return m_value[parNum];
  }

write()

void MuonResidualsFitter::write	(	FILE *	file,
		int	which = 0
	)

Definition at line 383 of file MuonResidualsFitter.cc.

References geometryDiff::file, mps_fire::i, ndata(), numResiduals(), residuals_begin(), residuals_end(), postprocess-scan-build::rows, and eostools::which().

Referenced by pkg.AbstractPkg::generate(), MuonResidualsTwoBin::write(), and querying.connection::write_and_commit().

                                                     {
  long rows = numResiduals();
  int cols = ndata();
  int whichcopy = which;

  fwrite(&rows, sizeof(long), 1, file);
  fwrite(&cols, sizeof(int), 1, file);
  fwrite(&whichcopy, sizeof(int), 1, file);

  double *likeAChecksum = new double[cols];
  double *likeAChecksum2 = new double[cols];
  for (int i = 0; i < cols; i++) {
    likeAChecksum[i] = 0.;
    likeAChecksum2[i] = 0.;
  }

  for (std::vector<double *>::const_iterator residual = residuals_begin(); residual != residuals_end(); ++residual) {
    fwrite((*residual), sizeof(double), cols, file);
    for (int i = 0; i < cols; i++) {
      if (fabs((*residual)[i]) > likeAChecksum[i])
        likeAChecksum[i] = fabs((*residual)[i]);
      if (fabs((*residual)[i]) < likeAChecksum2[i])
        likeAChecksum2[i] = fabs((*residual)[i]);
    }
  }  // end loop over residuals

  // the idea is that mal-formed doubles are likely to be huge values (or tiny values)
  // because the exponent gets screwed up; we want to check for that
  fwrite(likeAChecksum, sizeof(double), cols, file);
  fwrite(likeAChecksum2, sizeof(double), cols, file);

  delete[] likeAChecksum;
  delete[] likeAChecksum2;
}

Member Data Documentation

m_center

double MuonResidualsFitter::m_center[20]

protected

Definition at line 245 of file MuonResidualsFitter.h.

Referenced by selectPeakResiduals().

m_cov

TMatrixDSym MuonResidualsFitter::m_cov

protected

Definition at line 237 of file MuonResidualsFitter.h.

Referenced by covarianceElement(), covarianceMatrix(), and dofit().

m_error

std::vector<double> MuonResidualsFitter::m_error

protected

Definition at line 236 of file MuonResidualsFitter.h.

Referenced by dofit(), and errorerror().

m_fixed

std::vector<bool> MuonResidualsFitter::m_fixed

protected

Definition at line 228 of file MuonResidualsFitter.h.

Referenced by fix(), fixed(), and nfixed().

m_loglikelihood

double MuonResidualsFitter::m_loglikelihood

protected

Definition at line 238 of file MuonResidualsFitter.h.

Referenced by dofit(), and loglikelihood().

m_minHits

int MuonResidualsFitter::m_minHits

protected

Definition at line 225 of file MuonResidualsFitter.h.

Referenced by MuonResidualsAngleFitter::fit(), MuonResidualsPositionFitter::fit(), and MuonResidualsBfieldAngleFitter::fit().

m_parNum2InitValue

std::map<int, double> MuonResidualsFitter::m_parNum2InitValue

protected

Definition at line 229 of file MuonResidualsFitter.h.

Referenced by setInitialValue().

m_parNum2parIdx

std::map<int, int> MuonResidualsFitter::m_parNum2parIdx

protected

Definition at line 240 of file MuonResidualsFitter.h.

Referenced by parNum2parIdx().

m_printLevel

int MuonResidualsFitter::m_printLevel

protected

Definition at line 230 of file MuonResidualsFitter.h.

Referenced by setPrintLevel().

m_radii

double MuonResidualsFitter::m_radii[20]

protected

Definition at line 246 of file MuonResidualsFitter.h.

Referenced by selectPeakResiduals().

m_residuals

std::vector<double *> MuonResidualsFitter::m_residuals

protected

Definition at line 232 of file MuonResidualsFitter.h.

Referenced by computeHistogramRangeAndBinning(), correctBField(), eraseNotSelectedResiduals(), fill(), histogramChi2GaussianFit(), numResiduals(), read(), residuals_begin(), residuals_end(), selectPeakResiduals(), and selectPeakResiduals_simple().

m_residuals_ok

std::vector<bool> MuonResidualsFitter::m_residuals_ok

protected

Definition at line 233 of file MuonResidualsFitter.h.

Referenced by correctBField(), eraseNotSelectedResiduals(), fiducialCuts(), fill(), read(), selectedResidualsFlags(), selectPeakResiduals(), and selectPeakResiduals_simple().

m_residualsModel

int MuonResidualsFitter::m_residualsModel

protected

Definition at line 224 of file MuonResidualsFitter.h.

Referenced by MuonResidualsFitter(), and residualsModel().

m_strategy

int MuonResidualsFitter::m_strategy

protected

Definition at line 230 of file MuonResidualsFitter.h.

Referenced by dofit(), and setStrategy().

m_useResiduals

int MuonResidualsFitter::m_useResiduals

protected

Definition at line 226 of file MuonResidualsFitter.h.

Referenced by useRes().

m_value

std::vector<double> MuonResidualsFitter::m_value

protected

Definition at line 235 of file MuonResidualsFitter.h.

Referenced by dofit(), and value().

m_weightAlignment

bool MuonResidualsFitter::m_weightAlignment

protected

Definition at line 227 of file MuonResidualsFitter.h.

Referenced by MuonResiduals1DOFFitter::fit(), MuonResiduals1DOFFitter::plot(), MuonResiduals6DOFrphiFitter::plot(), MuonResiduals5DOFFitter::plot(), MuonResiduals6DOFFitter::plot(), MuonResiduals1DOFFitter::sumofweights(), MuonResiduals6DOFrphiFitter::sumofweights(), MuonResiduals5DOFFitter::sumofweights(), and MuonResiduals6DOFFitter::sumofweights().