test
CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
List of all members | Classes | Public Types | Public Member Functions | Protected Member Functions | Protected Attributes
MuonResidualsFitter Class Referenceabstract

#include <MuonResidualsFitter.h>

Inheritance diagram for MuonResidualsFitter:
MuonResiduals1DOFFitter MuonResiduals5DOFFitter MuonResiduals6DOFFitter MuonResiduals6DOFrphiFitter MuonResidualsAngleFitter MuonResidualsBfieldAngleFitter MuonResidualsPositionFitter

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
}
 

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 fill (double *residual)
 
virtual bool fit (Alignable *ali)=0
 
void fix (int parNum, bool val=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 setPrintLevel (int printLevel)
 
void setStrategy (int strategy)
 
virtual double sumofweights ()=0
 
virtual int type () const =0
 
int useRes () const
 
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, 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 80 of file MuonResidualsFitter.h.

Member Enumeration Documentation

anonymous enum
anonymous enum
anonymous enum

Constructor & Destructor Documentation

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

Definition at line 130 of file MuonResidualsFitter.h.

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

132  {
135  throw cms::Exception("MuonResidualsFitter") << "unrecognized residualsModel";
136  };
tuple weightAlignment
Definition: align_cfg.py:30
list useResiduals
Definition: align_cfg.py:36
virtual MuonResidualsFitter::~MuonResidualsFitter ( )
inlinevirtual

Definition at line 138 of file MuonResidualsFitter.h.

References residuals_begin(), and residuals_end().

139  {
140  for (std::vector<double*>::const_iterator residual = residuals_begin(); residual != residuals_end(); ++residual) {
141  delete [] (*residual);
142  }
143  }
std::vector< double * >::const_iterator residuals_end() const
std::vector< double * >::const_iterator residuals_begin() const

Member Function Documentation

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

Definition at line 455 of file MuonResidualsFitter.cc.

References gather_cfg::cout, data, i, m_residuals, gen::n, NULL, numResiduals(), funct::pow(), alignCSCRings::r, python.multivaluedict::sort(), and eostools::which().

Referenced by histogramChi2GaussianFit().

456 {
457  // first, make a numeric array while discarding some crazy outliers
458  double *data = new double[numResiduals()];
459  int n = 0;
460  for (std::vector<double*>::const_iterator r = m_residuals.begin(); r != m_residuals.end(); r++)
461  if (fabs((*r)[which])<50.)
462  {
463  data[n] = (*r)[which];
464  n++;
465  }
466 
467  // compute "3 normal sigma" and regular interquantile ranges
468  const int n_quantiles = 7;
469  double probabilities[n_quantiles] = {0.00135, 0.02275, 0.25, 0.5, 0.75, 0.97725, 0.99865}; // "3 normal sigma"
470  //double probabilities[n_quantiles] = {0.02275, 0.25, 0.75, 0.97725}; // "2 normal sigma"
471  double quantiles[n_quantiles];
472  std::sort(data, data + n);
473  TMath::Quantiles(n, n_quantiles, data, quantiles, probabilities, true, NULL, 7);
474  delete [] data;
475  double iqr = quantiles[4] - quantiles[2];
476 
477  // estimate optimal bin size according to Freedman-Diaconis rule
478  double hbin = 2 * iqr / pow( n, 1./3);
479 
480  a = quantiles[1];
481  b = quantiles[5];
482  nbins = (int) ( (b - a) / hbin + 3. ); // add extra safety margin of 3
483 
484  std::cout<<" quantiles: "; for (int i=0;i<n_quantiles;i++) std::cout<<quantiles[i]<<" "; std::cout<<std::endl;
485  //cout<<"n="<<select_count<<" quantiles ["<<quantiles[1]<<", "<<quantiles[2]<<"] IQR="<<iqr
486  // <<" full range=["<<minx<<","<<maxx<<"]"<<" 2 normal sigma quantile range = ["<<quantiles[0]<<", "<<quantiles[3]<<"]"<<endl;
487  std::cout<<" optimal h="<<hbin<<" nbins="<<nbins<<std::endl;
488 }
int i
Definition: DBlmapReader.cc:9
def which
Definition: eostools.py:335
#define NULL
Definition: scimark2.h:8
std::vector< double * > m_residuals
double b
Definition: hdecay.h:120
char data[epos_bytes_allocation]
Definition: EPOS_Wrapper.h:82
double a
Definition: hdecay.h:121
tuple cout
Definition: gather_cfg.py:121
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:40
virtual void MuonResidualsFitter::correctBField ( )
pure virtual
void MuonResidualsFitter::correctBField ( int  idx_momentum,
int  idx_q 
)
virtual

Definition at line 646 of file MuonResidualsFitter.cc.

References newFWLiteAna::bin, prof2calltree::count, gather_cfg::cout, alignCSCRings::e, i, j, m_residuals, m_residuals_ok, EnergyCorrector::pt, alignCSCRings::r, residuals_begin(), residuals_end(), findQualityFiles::size, and python.multivaluedict::sort().

647 {
648  const int Nbin = 17;
649  // find max 1/pt and bin width
650  double min_pt = 9999.;
651  for (std::vector<double*>::const_iterator r = residuals_begin(); r != residuals_end(); ++r)
652  {
653  double pt = fabs((*r)[idx_momentum]);
654  if (pt < min_pt) min_pt = pt;
655  }
656  min_pt -= 0.01; // to prevent bin # overflow
657  const double bin_width = 1./min_pt/Nbin;
658 
659  // fill indices of positive and negative charge residuals in each bin
660  std::vector<size_t> pos[Nbin], neg[Nbin], to_erase;
661  for (size_t i = 0; i < m_residuals.size(); i++)
662  {
663  if (!m_residuals_ok[i]) continue;
664  int bin = (int)floor(1./fabs(m_residuals[i][idx_momentum])/bin_width);
665  if (m_residuals[i][idx_q] > 0) pos[bin].push_back(i);
666  else neg[bin].push_back(i);
667  }
668 
669  // equalize pos and neg in each bin
670  for (int j = 0; j < Nbin; j++)
671  {
672  size_t psize = pos[j].size();
673  size_t nsize = neg[j].size();
674  if (psize == nsize) continue;
675 
676  std::set<int> idx_set; // use a set to collect certain number of unique random indices to erase
677  if (psize > nsize)
678  {
679  while (idx_set.size() < psize - nsize) idx_set.insert( gRandom->Integer(psize) );
680  for (std::set<int>::iterator it = idx_set.begin() ; it != idx_set.end(); it++ ) to_erase.push_back(pos[j][*it]);
681  }
682  else
683  {
684  while (idx_set.size() < nsize - psize) idx_set.insert( gRandom->Integer(nsize) );
685  for (std::set<int>::iterator it = idx_set.begin() ; it != idx_set.end(); it++ ) to_erase.push_back(neg[j][*it]);
686  }
687  }
688  // sort in descending order, so we safely go from higher to lower indices:
689  std::sort(to_erase.begin(), to_erase.end(), std::greater<int>() );
690  for (std::vector<size_t>::const_iterator e = to_erase.begin(); e != to_erase.end(); ++e)
691  {
692  m_residuals_ok[*e] = false;
693  //delete[] *(m_residuals.begin() + *e);
694  //m_residuals.erase(m_residuals.begin() + *e);
695  }
696 
697  std::vector<size_t> apos[Nbin], aneg[Nbin];
698  for (size_t i = 0; i < m_residuals.size(); i++)
699  {
700  if (!m_residuals_ok[i]) continue;
701  int bin = (int)floor(1./fabs(m_residuals[i][idx_momentum])/bin_width);
702  if (m_residuals[i][idx_q] > 0) apos[bin].push_back(i);
703  else aneg[bin].push_back(i);
704  }
705  for (int j = 0; j < Nbin; j++) std::cout << "bin " << j << ": [pos,neg] sizes before & after: ["
706  << pos[j].size() <<","<< neg[j].size()<<"] -> [" << apos[j].size() <<","<< aneg[j].size() << "]" << std::endl;
707  std::cout<<" N residuals "<<m_residuals.size()<<" -> "<<(size_t) std::count(m_residuals_ok.begin(), m_residuals_ok.end(), true)<<std::endl;
708 }
int i
Definition: DBlmapReader.cc:9
std::vector< double * > m_residuals
int j
Definition: DBlmapReader.cc:9
std::vector< bool > m_residuals_ok
std::vector< double * >::const_iterator residuals_end() const
tuple cout
Definition: gather_cfg.py:121
std::vector< double * >::const_iterator residuals_begin() const
tuple size
Write out results.
TMatrixDSym MuonResidualsFitter::correlationMatrix ( )
double MuonResidualsFitter::covarianceElement ( int  parNum1,
int  parNum2 
)
inline

Definition at line 191 of file MuonResidualsFitter.h.

References assert(), m_cov, npar(), and parNum2parIdx().

192  {
193  assert(0 <= parNum1 && parNum1 < npar());
194  assert(0 <= parNum2 && parNum2 < npar());
195  assert(m_cov.GetNcols() == npar()); // m_cov might have not yet been resized to account for proper #parameters
196  return m_cov(parNum2parIdx(parNum1), parNum2parIdx(parNum2));
197  }
int parNum2parIdx(int parNum)
virtual int npar()=0
assert(m_qm.get())
TMatrixDSym MuonResidualsFitter::covarianceMatrix ( )
inline

Definition at line 190 of file MuonResidualsFitter.h.

References m_cov.

190 {return m_cov;}
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 237 of file MuonResidualsFitter.cc.

References fitWZ::arglist, alignCSCRings::e, fcn(), fixed(), i, inform(), m_cov, m_error, m_loglikelihood, m_printLevel, m_strategy, m_value, MuonResidualsFitter_TMinuit, npar(), and findQualityFiles::v.

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

239 {
241 
242  MuonResidualsFitter_TMinuit = new TMinuit(npar());
244  MuonResidualsFitter_TMinuit->SetObjectFit(fitinfo);
247 
248  std::vector<int>::const_iterator iNum = parNum.begin();
249  std::vector<std::string>::const_iterator iName = parName.begin();
250  std::vector<double>::const_iterator istart = start.begin();
251  std::vector<double>::const_iterator istep = step.begin();
252  std::vector<double>::const_iterator ilow = low.begin();
253  std::vector<double>::const_iterator ihigh = high.begin();
254 
255  //MuonResidualsFitter_TMinuit->SetPrintLevel(-1);
256 
257  for (; iNum != parNum.end(); ++iNum, ++iName, ++istart, ++istep, ++ilow, ++ihigh)
258  {
259  MuonResidualsFitter_TMinuit->DefineParameter(*iNum, iName->c_str(), *istart, *istep, *ilow, *ihigh);
260  if (fixed(*iNum)) MuonResidualsFitter_TMinuit->FixParameter(*iNum);
261  }
262 
263  double arglist[10];
264  int ierflg;
265  int smierflg; //second MIGRAD ierflg
266 
267  // chi^2 errors should be 1.0, log-likelihood should be 0.5
268  for (int i = 0; i < 10; i++) arglist[i] = 0.;
269  arglist[0] = 0.5;
270  ierflg = 0;
271  smierflg = 0;
272  MuonResidualsFitter_TMinuit->mnexcm("SET ERR", arglist, 1, ierflg);
273  if (ierflg != 0) { delete MuonResidualsFitter_TMinuit; delete fitinfo; return false; }
274 
275  // set strategy = 2 (more refined fits)
276  for (int i = 0; i < 10; i++) arglist[i] = 0.;
277  arglist[0] = m_strategy;
278  ierflg = 0;
279  MuonResidualsFitter_TMinuit->mnexcm("SET STR", arglist, 1, ierflg);
280  if (ierflg != 0) { delete MuonResidualsFitter_TMinuit; delete fitinfo; return false; }
281 
282  bool try_again = false;
283 
284  // minimize
285  for (int i = 0; i < 10; i++) arglist[i] = 0.;
286  arglist[0] = 50000;
287  ierflg = 0;
288  MuonResidualsFitter_TMinuit->mnexcm("MIGRAD", arglist, 1, ierflg);
289  if (ierflg != 0) try_again = true;
290 
291  // just once more, if needed (using the final Minuit parameters from the failed fit; often works)
292  if (try_again)
293  {
294  for (int i = 0; i < 10; i++) arglist[i] = 0.;
295  arglist[0] = 50000;
296  MuonResidualsFitter_TMinuit->mnexcm("MIGRAD", arglist, 1, smierflg);
297  }
298 
299  Double_t fmin, fedm, errdef;
300  Int_t npari, nparx, istat;
301  MuonResidualsFitter_TMinuit->mnstat(fmin, fedm, errdef, npari, nparx, istat);
302 
303  if (istat != 3)
304  {
305  for (int i = 0; i < 10; i++) arglist[i] = 0.;
306  ierflg = 0;
307  MuonResidualsFitter_TMinuit->mnexcm("HESSE", arglist, 0, ierflg);
308  }
309 
310  // read-out the results
311  m_loglikelihood = -fmin;
312 
313  m_value.clear();
314  m_error.clear();
315  for (int i = 0; i < npar(); i++)
316  {
317  double v, e;
318  MuonResidualsFitter_TMinuit->GetParameter(i, v, e);
319  m_value.push_back(v);
320  m_error.push_back(e);
321  }
322  m_cov.ResizeTo(npar(),npar());
323  MuonResidualsFitter_TMinuit->mnemat( m_cov.GetMatrixArray(), npar());
324 
326  delete fitinfo;
327  if (smierflg != 0) return false;
328  return true;
329 }
tuple arglist
Definition: fitWZ.py:38
int i
Definition: DBlmapReader.cc:9
tuple start
Check for commandline option errors.
Definition: dqm_diff.py:58
virtual void inform(TMinuit *tMinuit)=0
virtual int npar()=0
std::vector< double > m_error
std::vector< double > m_value
bool fixed(int parNum)
static TMinuit * MuonResidualsFitter_TMinuit
void fcn(int &, double *, double &, double *, int)
void MuonResidualsFitter::eraseNotSelectedResiduals ( )

Definition at line 711 of file MuonResidualsFitter.cc.

References prof2calltree::count, gather_cfg::cout, i, m_residuals, m_residuals_ok, and tmp.

Referenced by MuonResidualsTwoBin::eraseNotSelectedResiduals().

712 {
713  // it should probably be faster then doing erase
714  size_t n_ok = (size_t) std::count(m_residuals_ok.begin(), m_residuals_ok.end(), true);
715  std::vector<double*> tmp(n_ok, 0);
716  std::cout << "residuals sizes: all=" << m_residuals.size()<<" good="<<n_ok<<std::endl;
717  int iok=0;
718  for (size_t i = 0; i < m_residuals.size(); i++)
719  {
720  if (!m_residuals_ok[i]) continue;
721  tmp[iok++] = m_residuals[i];
722  }
723  m_residuals.swap(tmp);
724 
725  std::vector<bool> tmp_ok(n_ok, true);
726  m_residuals_ok.swap(tmp_ok);
727 
728  std::cout << "residuals size after eraseNotSelectedResiduals =" << m_residuals.size()<<" ok size="<<m_residuals_ok.size()<<std::endl;
729 }
int i
Definition: DBlmapReader.cc:9
std::vector< double * > m_residuals
std::vector< bool > m_residuals_ok
std::vector< std::vector< double > > tmp
Definition: MVATrainer.cc:100
tuple cout
Definition: gather_cfg.py:121
double MuonResidualsFitter::errorerror ( int  parNum)
inline

Definition at line 184 of file MuonResidualsFitter.h.

References assert(), m_error, and npar().

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

184 { assert(0 <= parNum && parNum < npar()); return m_error[parNum]; }
virtual int npar()=0
std::vector< double > m_error
assert(m_qm.get())
void MuonResidualsFitter::fill ( double *  residual)
inline

Definition at line 173 of file MuonResidualsFitter.h.

References m_residuals, and m_residuals_ok.

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

174  {
175  m_residuals.push_back(residual);
176  m_residuals_ok.push_back(true);
177  }
std::vector< double * > m_residuals
std::vector< bool > m_residuals_ok
virtual bool MuonResidualsFitter::fit ( Alignable ali)
pure virtual
void MuonResidualsFitter::fix ( int  parNum,
bool  val = true 
)
inline

Definition at line 153 of file MuonResidualsFitter.h.

References assert(), m_fixed, and npar().

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

154  {
155  assert(0 <= parNum && parNum < npar());
156  if (m_fixed.size() == 0) m_fixed.resize(npar(), false);
157  m_fixed[parNum] = val;
158  }
virtual int npar()=0
assert(m_qm.get())
std::vector< bool > m_fixed
bool MuonResidualsFitter::fixed ( int  parNum)
inline

Definition at line 160 of file MuonResidualsFitter.h.

References assert(), m_fixed, and npar().

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

161  {
162  assert(0 <= parNum && parNum < npar());
163  if (m_fixed.size() == 0) return false;
164  else return m_fixed[parNum];
165  }
virtual int npar()=0
assert(m_qm.get())
std::vector< bool > m_fixed
void MuonResidualsFitter::histogramChi2GaussianFit ( int  which,
double &  fit_mean,
double &  fit_sigma 
)

Definition at line 491 of file MuonResidualsFitter.cc.

References a, b, computeHistogramRangeAndBinning(), gather_cfg::cout, python.connectstrParser::f1, estimatePileup::hist, m_residuals, pileupCalc::nbins, and alignCSCRings::r.

Referenced by selectPeakResiduals_simple().

492 {
493  int nbins;
494  double a, b;
495  computeHistogramRangeAndBinning(which, nbins, a, b);
496  if (a==b || a > b) { fit_mean = a; fit_sigma = 0; return; }
497 
498  TH1D *hist = new TH1D("htmp", "", nbins, a, b);
499  for (std::vector<double*>::const_iterator r = m_residuals.begin(); r != m_residuals.end(); ++r) hist->Fill( (*r)[which] );
500 
501  // do simple chi2 gaussian fit
502  TF1 *f1= new TF1("f1","gaus", a, b);
503  f1->SetParameter(0, hist->GetEntries());
504  f1->SetParameter(1, 0);
505  f1->SetParameter(2, hist->GetRMS());
506  hist->Fit(f1,"RQ");
507 
508  fit_mean = f1->GetParameter(1);
509  fit_sigma = f1->GetParameter(2);
510  std::cout<<" h("<<nbins<<","<<a<<","<<b<<") mu="<<fit_mean<<" sig="<<fit_sigma<<std::endl;
511 
512  delete f1;
513  delete hist;
514 }
def which
Definition: eostools.py:335
void computeHistogramRangeAndBinning(int which, int &nbins, double &a, double &b)
std::vector< double * > m_residuals
double b
Definition: hdecay.h:120
double a
Definition: hdecay.h:121
tuple cout
Definition: gather_cfg.py:121
virtual void MuonResidualsFitter::inform ( TMinuit *  tMinuit)
protectedpure virtual
void MuonResidualsFitter::initialize_table ( )
protected

Definition at line 166 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, and residualsModel().

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

167 {
170 
171  std::ifstream convolution_table("convolution_table.txt");
172  if (convolution_table.is_open())
173  {
174  int numgsbins = 0;
175  int numtsbins = 0;
176  double tsbinsize = 0.;
177  double gsbinsize = 0.;
178 
179  convolution_table >> numgsbins >> numtsbins >> tsbinsize >> gsbinsize;
180  if (numgsbins != MuonResidualsFitter_numgsbins || numtsbins != MuonResidualsFitter_numtsbins ||
182  {
183  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";
184  }
185 
186  for (int gsbin = 0; gsbin < MuonResidualsFitter_numgsbins; gsbin++)
187  {
188  for (int tsbin = 0; tsbin < MuonResidualsFitter_numtsbins; tsbin++)
189  {
190  int read_gsbin = 0;
191  int read_tsbin = 0;
192  double val = 0.;
193 
194  convolution_table >> read_gsbin >> read_tsbin >> val;
195  if (read_gsbin != gsbin || read_tsbin != tsbin)
196  {
197  throw cms::Exception("MuonResidualsFitter") << "convolution_table.txt is out of order. Throw it away and let the fitter re-create the file.\n";
198  }
199  MuonResidualsFitter_lookup_table[gsbin][tsbin] = val;
200  }
201  }
202  convolution_table.close();
203  }
204  else
205  {
206  std::ofstream convolution_table2("convolution_table.txt");
207 
208  if (!convolution_table2.is_open()) throw cms::Exception("MuonResidualsFitter") << "Couldn't write to file convolution_table.txt\n";
209 
210  convolution_table2 << MuonResidualsFitter_numgsbins << " " << MuonResidualsFitter_numtsbins << " " << MuonResidualsFitter_tsbinsize << " " << MuonResidualsFitter_gsbinsize << std::endl;
211 
212  std::cout << "Initializing convolution look-up table (takes a few minutes)..." << std::endl;
213 
214  for (int gsbin = 0; gsbin < MuonResidualsFitter_numgsbins; gsbin++)
215  {
216  double gammaoversigma = double(gsbin) * MuonResidualsFitter_gsbinsize;
217  std::cout << " gsbin " << gsbin << "/" << MuonResidualsFitter_numgsbins << std::endl;
218  for (int tsbin = 0; tsbin < MuonResidualsFitter_numtsbins; tsbin++)
219  {
220  double toversigma = double(tsbin) * MuonResidualsFitter_tsbinsize;
221 
222  // 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)
223  MuonResidualsFitter_lookup_table[gsbin][tsbin] = MuonResidualsFitter_compute_log_convolution(toversigma, gammaoversigma);
224 
225  // <10% errors with max=20, step=0.005, power=4 (faster computation for testing)
226  // MuonResidualsFitter_lookup_table[gsbin][tsbin] = MuonResidualsFitter_compute_log_convolution(toversigma, gammaoversigma, 100., 0.005, 4.);
227 
228  convolution_table2 << gsbin << " " << tsbin << " " << MuonResidualsFitter_lookup_table[gsbin][tsbin] << std::endl;
229  }
230  }
231  convolution_table2.close();
232  std::cout << "Initialization done!" << std::endl;
233  }
234 }
const int MuonResidualsFitter_numgsbins
double MuonResidualsFitter_compute_log_convolution(double toversigma, double gammaoversigma, double max=1000., double step=0.001, double power=4.)
const int MuonResidualsFitter_numtsbins
const double MuonResidualsFitter_gsbinsize
const double MuonResidualsFitter_tsbinsize
bool MuonResidualsFitter_table_initialized
double MuonResidualsFitter_lookup_table[MuonResidualsFitter_numgsbins][MuonResidualsFitter_numtsbins]
tuple cout
Definition: gather_cfg.py:121
double MuonResidualsFitter::loglikelihood ( )
inline

Definition at line 200 of file MuonResidualsFitter.h.

References m_loglikelihood.

Referenced by MuonResidualsTwoBin::loglikelihood().

200 { return m_loglikelihood; }
virtual int MuonResidualsFitter::ndata ( )
pure virtual
int MuonResidualsFitter::nfixed ( )
inline

Definition at line 166 of file MuonResidualsFitter.h.

References prof2calltree::count, and m_fixed.

166 { return std::count(m_fixed.begin(), m_fixed.end(), true); }
std::vector< bool > m_fixed
virtual int MuonResidualsFitter::npar ( )
pure virtual
long MuonResidualsFitter::numResiduals ( ) const
inline
long MuonResidualsFitter::numsegments ( )
inline

Definition at line 202 of file MuonResidualsFitter.h.

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

Referenced by MuonResidualsTwoBin::numsegments().

203  {
204  long num = 0;
205  for (std::vector<double*>::const_iterator resiter = residuals_begin(); resiter != residuals_end(); ++resiter) num++;
206  return num;
207  }
std::vector< double * >::const_iterator residuals_end() const
std::vector< double * >::const_iterator residuals_begin() const
int MuonResidualsFitter::parNum2parIdx ( int  parNum)
inline

Definition at line 188 of file MuonResidualsFitter.h.

References m_parNum2parIdx.

Referenced by covarianceElement().

188 { return m_parNum2parIdx[parNum];}
std::map< int, int > m_parNum2parIdx
virtual double MuonResidualsFitter::plot ( std::string  name,
TFileDirectory dir,
Alignable ali 
)
pure virtual
void MuonResidualsFitter::plotsimple ( std::string  name,
TFileDirectory dir,
int  which,
double  multiplier 
)

Definition at line 425 of file MuonResidualsFitter.cc.

References estimatePileup::hist, TFileDirectory::make(), alignCSCRings::r, residuals_begin(), residuals_end(), and svgfig::window().

Referenced by MuonResidualsTwoBin::plotsimple().

426 {
427  double window = 100.;
428  if (which == 0) window = 2.*30.;
429  else if (which == 1) window = 2.*30.;
430  else if (which == 2) window = 2.*20.;
431  else if (which == 3) window = 2.*50.;
432 
433  TH1F *hist = dir->make<TH1F>(name.c_str(), "", 200, -window, window);
434  for (std::vector<double*>::const_iterator r = residuals_begin(); r != residuals_end(); ++r) hist->Fill(multiplier * (*r)[which]);
435 }
def window
Definition: svgfig.py:642
def which
Definition: eostools.py:335
T * make(const Args &...args) const
make new ROOT object
std::vector< double * >::const_iterator residuals_end() const
std::vector< double * >::const_iterator residuals_begin() const
void MuonResidualsFitter::plotweighted ( std::string  name,
TFileDirectory dir,
int  which,
int  whichredchi2,
double  multiplier 
)

Definition at line 438 of file MuonResidualsFitter.cc.

References estimatePileup::hist, TFileDirectory::make(), alignCSCRings::r, residuals_begin(), residuals_end(), puppiForMET_cff::weight, and svgfig::window().

Referenced by MuonResidualsTwoBin::plotweighted().

439 {
440  double window = 100.;
441  if (which == 0) window = 2.*30.;
442  else if (which == 1) window = 2.*30.;
443  else if (which == 2) window = 2.*20.;
444  else if (which == 3) window = 2.*50.;
445 
446  TH1F *hist = dir->make<TH1F>(name.c_str(), "", 200, -window, window);
447  for (std::vector<double*>::const_iterator r = residuals_begin(); r != residuals_end(); ++r)
448  {
449  double weight = 1./(*r)[whichredchi2];
450  if (TMath::Prob(1./weight*12, 12) < 0.99) hist->Fill(multiplier * (*r)[which], weight);
451  }
452 }
def window
Definition: svgfig.py:642
def which
Definition: eostools.py:335
T * make(const Args &...args) const
make new ROOT object
std::vector< double * >::const_iterator residuals_end() const
std::vector< double * >::const_iterator residuals_begin() const
void MuonResidualsFitter::read ( FILE *  file,
int  which = 0 
)

Definition at line 369 of file MuonResidualsFitter.cc.

References alignCSCRings::e, Exception, fill(), i, m_residuals, m_residuals_ok, ndata(), numResiduals(), and tablePrinter::rows.

Referenced by edmIntegrityCheck.PublishToFileSystem::get(), Vispa.Plugins.EdmBrowser.EdmDataAccessor.EdmDataAccessor::goto(), MuonResidualsTwoBin::read(), and Vispa.Plugins.EdmBrowser.EdmDataAccessor.EdmDataAccessor::setFilterBranches().

370 {
371  long rows = -100;
372  int cols = -100;
373  int readwhich = -100;
374 
375  fread(&rows, sizeof(long), 1, file);
376  fread(&cols, sizeof(int), 1, file);
377  fread(&readwhich, sizeof(int), 1, file);
378 
379  if (cols != ndata() || rows < 0 || readwhich != which)
380  {
381  throw cms::Exception("MuonResidualsFitter") << "temporary file is corrupted (which = " << which << " readwhich = " << readwhich << " rows = " << rows << " cols = " << cols << ")\n";
382  }
383 
384  double *likeAChecksum = new double[cols];
385  double *likeAChecksum2 = new double[cols];
386  for (int i = 0; i < cols; i++)
387  {
388  likeAChecksum[i] = 0.;
389  likeAChecksum2[i] = 0.;
390  }
391 
392  m_residuals.reserve(rows);
393  for (long row = 0; row < rows; row++)
394  {
395  double *residual = new double[cols];
396  fread(residual, sizeof(double), cols, file);
397  fill(residual);
398  for (int i = 0; i < cols; i++)
399  {
400  if (fabs(residual[i]) > likeAChecksum[i]) likeAChecksum[i] = fabs(residual[i]);
401  if (fabs(residual[i]) < likeAChecksum2[i]) likeAChecksum2[i] = fabs(residual[i]);
402  }
403  } // end loop over records in file
404 
405  double *readChecksum = new double[cols];
406  double *readChecksum2 = new double[cols];
407  fread(readChecksum, sizeof(double), cols, file);
408  fread(readChecksum2, sizeof(double), cols, file);
409 
410  for (int i = 0; i < cols; i++)
411  {
412  if (fabs(likeAChecksum[i] - readChecksum[i]) > 1e-10 || fabs(1./likeAChecksum2[i] - 1./readChecksum2[i]) > 1e10)
413  {
414  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";
415  }
416  }
417 
418  m_residuals_ok.resize(numResiduals(), true);
419 
420  delete [] likeAChecksum;
421  delete [] likeAChecksum2;
422 }
int i
Definition: DBlmapReader.cc:9
def which
Definition: eostools.py:335
void fill(double *residual)
std::vector< double * > m_residuals
std::vector< bool > m_residuals_ok
virtual int ndata()=0
std::vector<double*>::const_iterator MuonResidualsFitter::residuals_begin ( ) const
inline
std::vector<double*>::const_iterator MuonResidualsFitter::residuals_end ( ) const
inline
int MuonResidualsFitter::residualsModel ( ) const
inline
std::vector<bool>& MuonResidualsFitter::selectedResidualsFlags ( )
inline
void MuonResidualsFitter::selectPeakResiduals ( double  nsigma,
int  nvar,
int *  vars 
)

Definition at line 559 of file MuonResidualsFitter.cc.

References assert(), prof2calltree::count, gather_cfg::cout, data, HLT_25ns14e33_v1_cff::distance, i, m_residuals, m_residuals_ok, numResiduals(), alignCSCRings::r, mathSSE::sqrt(), and findQualityFiles::v.

Referenced by MuonResidualsTwoBin::selectPeakResiduals().

560 {
561  //std::cout<<"doing selectpeakresiduals: nsig="<<nsigma<<" nvar="<<nvar<<" vars=";
562  for (int i=0; i<nvar; ++i) std::cout<<vars[i]<<" ";
563  std::cout<<std::endl;
564 
565  // does not make sense for small statistics
566  if (numResiduals()<50) return;
567 
568  size_t nbefore = numResiduals();
569  std::cout<<" N residuals "<<nbefore<<" ~ "<<(size_t) std::count(m_residuals_ok.begin(), m_residuals_ok.end(), true)<<std::endl;
570  //just to be sure (can't see why it might ever be more then 10)
571  assert(nvar<=10 && nvar>0);
572 
573  std::vector<double*>::iterator r = m_residuals.begin();
574 
575  // it's awkward, but the 1D case has to be handled separately
576  if (nvar==1)
577  {
578  // get robust estimates for the peak and sigma
579  double *data = new double[nbefore];
580  for (size_t i = 0; i < nbefore; i++) data[i] = m_residuals[i][ vars[0] ];
581  double peak, sigma;
582  TRobustEstimator re;
583  re.EvaluateUni(nbefore, data, peak, sigma);
584 
585  // filter out residuals that are more then nsigma away from the peak
586  while (r != m_residuals.end())
587  {
588  double distance = fabs( ((*r)[ vars[0] ] - peak)/sigma );
589  if (distance <= nsigma) ++r;
590  else
591  {
592  m_residuals_ok[r - m_residuals.begin()] = false;
593  ++r;
594  //delete [] (*r);
595  //r = m_residuals.erase(r);
596  }
597  }
598  std::cout<<" N residuals "<<nbefore<<" -> "<<numResiduals()<<std::endl;
599  return;
600  } // end 1D case
601 
602  // initialize and run the robust estimator for D>1
603  TRobustEstimator re(nbefore, nvar);
604  r = m_residuals.begin();
605  while (r != m_residuals.end())
606  {
607  double *row = new double[nvar];
608  for (int v = 0; v<nvar; v++) row[v] = (*r)[ vars[v] ];
609  re.AddRow(row);
610  delete[] row;
611  ++r;
612  }
613  re.Evaluate();
614 
615  // get nvar-dimensional ellipsoid center & covariance
616  TVectorD M(nvar);
617  re.GetMean(M);
618  TMatrixDSym Cov(nvar);
619  re.GetCovariance(Cov);
620  Cov.Invert();
621 
622  // calculate the normalized radius for this nvar-dimensional ellipsoid from a 1D-Gaussian nsigma equivalent distance
623  double conf_1d = TMath::Erf(nsigma/sqrt(2));
624  double surf_radius = sqrt(TMath::ChisquareQuantile(conf_1d, nvar));
625 
626  // filter out residuals that are outside of the covariance ellipsoid with the above normalized radius
627  r = m_residuals.begin();
628  while (r != m_residuals.end())
629  {
630  TVectorD res(nvar);
631  for (int v = 0; v<nvar; v++) res[v] = (*r)[ vars[v] ];
632  double distance = sqrt( Cov.Similarity(res - M) );
633  if (distance <= surf_radius) ++r;
634  else
635  {
636  m_residuals_ok[r - m_residuals.begin()] = false;
637  ++r;
638  //delete [] (*r);
639  //r = m_residuals.erase(r);
640  }
641  }
642  std::cout<<" N residuals "<<nbefore<<" -> "<<(size_t) std::count(m_residuals_ok.begin(), m_residuals_ok.end(), true)<<std::endl;
643 }
int i
Definition: DBlmapReader.cc:9
assert(m_qm.get())
std::vector< double * > m_residuals
T sqrt(T t)
Definition: SSEVec.h:48
std::vector< bool > m_residuals_ok
char data[epos_bytes_allocation]
Definition: EPOS_Wrapper.h:82
tuple cout
Definition: gather_cfg.py:121
void MuonResidualsFitter::selectPeakResiduals_simple ( double  nsigma,
int  nvar,
int *  vars 
)

Definition at line 518 of file MuonResidualsFitter.cc.

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

519 {
520  // does not make sense for small statistics
521  if (numResiduals()<25) return;
522 
523  int nbefore = numResiduals();
524 
525  //just to be sure (can't see why it might ever be more then 10)
526  assert(nvar<=10);
527 
528  // estimate nvar-D ellipsoid center & axes
529  for (int v = 0; v<nvar; v++)
530  {
531  int which = vars[v];
532  histogramChi2GaussianFit(which, m_center[which], m_radii[which]);
533  m_radii[which] = nsigma * m_radii[which];
534  }
535 
536  // filter out residuals that don't fit into the ellipsoid
537  std::vector<double*>::iterator r = m_residuals.begin();
538  while (r != m_residuals.end())
539  {
540  double ellipsoid_sum = 0;
541  for (int v = 0; v<nvar; v++)
542  {
543  int which = vars[v];
544  if (m_radii[which] == 0.) continue;
545  ellipsoid_sum += pow( ( (*r)[which] - m_center[which]) / m_radii[which] , 2);
546  }
547  if (ellipsoid_sum <= 1.) ++r;
548  else
549  {
550  delete [] (*r);
551  r = m_residuals.erase(r);
552  }
553  }
554  std::cout<<" N residuals "<<nbefore<<" -> "<<numResiduals()<<std::endl;
555 }
void histogramChi2GaussianFit(int which, double &fit_mean, double &fit_sigma)
def which
Definition: eostools.py:335
assert(m_qm.get())
std::vector< double * > m_residuals
tuple cout
Definition: gather_cfg.py:121
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:40
void MuonResidualsFitter::setPrintLevel ( int  printLevel)
inline

Definition at line 168 of file MuonResidualsFitter.h.

References m_printLevel.

Referenced by MuonResidualsTwoBin::setPrintLevel().

168 { m_printLevel = printLevel; }
void MuonResidualsFitter::setStrategy ( int  strategy)
inline

Definition at line 169 of file MuonResidualsFitter.h.

References m_strategy.

Referenced by MuonResidualsTwoBin::setStrategy().

169 { m_strategy = strategy; }
virtual double MuonResidualsFitter::sumofweights ( )
pure virtual
virtual int MuonResidualsFitter::type ( ) const
pure virtual
int MuonResidualsFitter::useRes ( ) const
inline
double MuonResidualsFitter::value ( int  parNum)
inline

Definition at line 183 of file MuonResidualsFitter.h.

References 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(), Vispa.Views.PropertyView.FileProperty::buttonClicked(), Types.string::configValue(), Types.FileInPath::configValue(), Mixins.UsingBlock::dumpPython(), Types.int32::insertInto(), Types.uint32::insertInto(), Types.int64::insertInto(), Types.uint64::insertInto(), Mixins.UsingBlock::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(), Vispa.Views.PropertyView.FileProperty::labelDoubleClicked(), MuonResiduals1DOFFitter::plot(), MuonResidualsAngleFitter::plot(), MuonResidualsBfieldAngleFitter::plot(), MuonResidualsPositionFitter::plot(), MuonResiduals5DOFFitter::plot(), MuonResiduals6DOFFitter::plot(), MuonResiduals6DOFrphiFitter::plot(), and MuonResidualsTwoBin::value().

183 { assert(0 <= parNum && parNum < npar()); return m_value[parNum]; }
virtual int npar()=0
assert(m_qm.get())
std::vector< double > m_value
void MuonResidualsFitter::write ( FILE *  file,
int  which = 0 
)

Definition at line 332 of file MuonResidualsFitter.cc.

References i, ndata(), numResiduals(), residuals_begin(), residuals_end(), tablePrinter::rows, and eostools::which().

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

333 {
334  long rows = numResiduals();
335  int cols = ndata();
336  int whichcopy = which;
337 
338  fwrite(&rows, sizeof(long), 1, file);
339  fwrite(&cols, sizeof(int), 1, file);
340  fwrite(&whichcopy, sizeof(int), 1, file);
341 
342  double *likeAChecksum = new double[cols];
343  double *likeAChecksum2 = new double[cols];
344  for (int i = 0; i < cols; i++)
345  {
346  likeAChecksum[i] = 0.;
347  likeAChecksum2[i] = 0.;
348  }
349 
350  for (std::vector<double*>::const_iterator residual = residuals_begin(); residual != residuals_end(); ++residual)
351  {
352  fwrite((*residual), sizeof(double), cols, file);
353  for (int i = 0; i < cols; i++)
354  {
355  if (fabs((*residual)[i]) > likeAChecksum[i]) likeAChecksum[i] = fabs((*residual)[i]);
356  if (fabs((*residual)[i]) < likeAChecksum2[i]) likeAChecksum2[i] = fabs((*residual)[i]);
357  }
358  } // end loop over residuals
359 
360  // the idea is that mal-formed doubles are likely to be huge values (or tiny values)
361  // because the exponent gets screwed up; we want to check for that
362  fwrite(likeAChecksum, sizeof(double), cols, file);
363  fwrite(likeAChecksum2, sizeof(double), cols, file);
364 
365  delete [] likeAChecksum;
366  delete [] likeAChecksum2;
367 }
int i
Definition: DBlmapReader.cc:9
def which
Definition: eostools.py:335
virtual int ndata()=0
std::vector< double * >::const_iterator residuals_end() const
std::vector< double * >::const_iterator residuals_begin() const

Member Data Documentation

double MuonResidualsFitter::m_center[20]
protected

Definition at line 271 of file MuonResidualsFitter.h.

Referenced by selectPeakResiduals_simple().

TMatrixDSym MuonResidualsFitter::m_cov
protected

Definition at line 263 of file MuonResidualsFitter.h.

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

std::vector<double> MuonResidualsFitter::m_error
protected

Definition at line 262 of file MuonResidualsFitter.h.

Referenced by dofit(), and errorerror().

std::vector<bool> MuonResidualsFitter::m_fixed
protected

Definition at line 255 of file MuonResidualsFitter.h.

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

double MuonResidualsFitter::m_loglikelihood
protected

Definition at line 264 of file MuonResidualsFitter.h.

Referenced by dofit(), and loglikelihood().

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

Definition at line 266 of file MuonResidualsFitter.h.

Referenced by parNum2parIdx().

int MuonResidualsFitter::m_printLevel
protected

Definition at line 256 of file MuonResidualsFitter.h.

Referenced by dofit(), and setPrintLevel().

double MuonResidualsFitter::m_radii[20]
protected

Definition at line 272 of file MuonResidualsFitter.h.

Referenced by selectPeakResiduals_simple().

std::vector<double*> MuonResidualsFitter::m_residuals
protected
std::vector<bool> MuonResidualsFitter::m_residuals_ok
protected
int MuonResidualsFitter::m_residualsModel
protected

Definition at line 251 of file MuonResidualsFitter.h.

Referenced by MuonResidualsFitter(), and residualsModel().

int MuonResidualsFitter::m_strategy
protected

Definition at line 256 of file MuonResidualsFitter.h.

Referenced by dofit(), and setStrategy().

int MuonResidualsFitter::m_useResiduals
protected

Definition at line 253 of file MuonResidualsFitter.h.

Referenced by useRes().

std::vector<double> MuonResidualsFitter::m_value
protected

Definition at line 261 of file MuonResidualsFitter.h.

Referenced by dofit(), and value().

bool MuonResidualsFitter::m_weightAlignment
protected