MuonResidualsFitter.h

Go to the documentation of this file.

#ifndef Alignment_MuonAlignmentAlgorithms_MuonResidualsFitter_H
#define Alignment_MuonAlignmentAlgorithms_MuonResidualsFitter_H

#ifndef STANDALONE_FITTER
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "CommonTools/UtilAlgos/interface/TFileService.h"
#include "Alignment/CommonAlignment/interface/Alignable.h"
#endif

#include "TMinuit.h"
#include "TH1F.h"
#include "TProfile.h"
#include "TF1.h"
#include "TMath.h"
#include "TRandom3.h"
#include "TMatrixDSym.h"

#include <cstdio>
#include <iostream>
#include <string>
#include <sstream>
#include <map>

// some mock classes for the case if we want to compile fitters with pure root outside of CMSSW
#ifdef STANDALONE_FITTER

#include <cassert>

class Alignable {
public:
  struct Surface {
    double width() { return 300.; }
    double length() { return 300.; }
  };
  Surface s;
  Surface &surface() { return s; }
};

class TFileDirectory {
public:
  template <typename T, typename A1, typename A2, typename A3, typename A4, typename A5>
  T *make(const A1 &a1, const A2 &a2, const A3 &a3, const A4 &a4, const A5 &a5) const {
    T *t = new T(a1, a2, a3, a4, a5);
    return t;
  }
  template <typename T, typename A1, typename A2, typename A3, typename A4, typename A5, typename A6, typename A7>
  T *make(const A1 &a1, const A2 &a2, const A3 &a3, const A4 &a4, const A5 &a5, const A6 &a6, const A7 &a7) const {
    T *t = new T(a1, a2, a3, a4, a5, a6, a7);
    return t;
  }
  template <typename T, typename A1, typename A2, typename A3, typename A4, typename A5, typename A6, typename A7, typename A8>
  T *make(const A1 &a1, const A2 &a2, const A3 &a3, const A4 &a4, const A5 &a5, const A6 &a6, const A7 &a7, const A8 &a8)
      const {
    T *t = new T(a1, a2, a3, a4, a5, a6, a7, a8);
    return t;
  }
};

#include <exception>
namespace cms {
  class Exception : public std::exception {
  public:
    Exception(const std::string &s) : name(s) {}
    ~Exception() throw() {}
    std::string name;
    template <class T>
    Exception &operator<<(const T &what) {
      std::cout << name << " exception: " << what << std::endl;
      return *this;
    }
  };
}  // namespace cms

#endif  // STANDALONE_FITTER

class MuonResidualsFitter {
public:
  enum { kPureGaussian, kPowerLawTails, kROOTVoigt, kGaussPowerTails, kPureGaussian2D };

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

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

  struct MuonAlignmentTreeRow {
    Bool_t is_plus;   // for +- endcap
    Bool_t is_dt;     // DT or CSC
    UChar_t station;  // 8bit uint
    Char_t ring_wheel;
    UChar_t sector;
    Float_t res_x;
    Float_t res_y;
    Float_t res_slope_x;
    Float_t res_slope_y;
    Float_t pos_x;
    Float_t pos_y;
    Float_t angle_x;
    Float_t angle_y;
    Float_t pz;
    Float_t pt;
    Char_t q;
    Bool_t select;
  };

  MuonResidualsFitter(int residualsModel, int minHits, int useResiduals, bool weightAlignment = true);
  virtual ~MuonResidualsFitter();

  virtual int type() const = 0;
  virtual int npar() = 0;
  virtual int ndata() = 0;

  int useRes(int pattern = -1) {
    if (pattern >= 0)
      m_useResiduals = pattern;
    return m_useResiduals;
  }
  int residualsModel() const { return m_residualsModel; }
  long numResiduals() const { return m_residuals.size(); }

  void fix(int parNum, bool dofix = true);
  bool fixed(int parNum);

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

  void setPrintLevel(int printLevel) { m_printLevel = printLevel; }
  void setStrategy(int strategy) { m_strategy = strategy; }

  void setInitialValue(int parNum, double value) { m_parNum2InitValue[parNum] = value; }

  // an array of the actual residual and associated baggage (qoverpt, trackangle, trackposition)
  // arrays passed to fill() are "owned" by MuonResidualsFitter: MuonResidualsFitter will delete them, don't do it yourself!
  void fill(double *residual);

  // this block of results is only valid if fit() returns true
  // also gamma is only valid if the model is kPowerLawTails or kROOTVoigt
  virtual bool fit(Alignable *ali) = 0;

  double value(int parNum) {
    assert(0 <= parNum && parNum < npar());
    return m_value[parNum];
  }
  double errorerror(int parNum) {
    assert(0 <= parNum && parNum < npar());
    return m_error[parNum];
  }

  // parNum corresponds to the parameter number defined by enums in specific fitters
  // parIdx is a continuous index in a set of parameters
  int parNum2parIdx(int parNum) { return m_parNum2parIdx[parNum]; }

  TMatrixDSym covarianceMatrix() { return m_cov; }
  double covarianceElement(int parNum1, int parNum2);
  TMatrixDSym correlationMatrix();

  double loglikelihood() { return m_loglikelihood; }

  long numsegments();

  virtual double sumofweights() = 0;

  // demonstration plots; return reduced chi**2
  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);

#ifdef STANDALONE_FITTER
  Alignable m_ali;
  TFileDirectory m_dir;
  bool fit() { return fit(&m_ali); }
  virtual double plot(std::string &name) { return plot(name, &m_dir, &m_ali); }
  void plotsimple(std::string &name, int which, double multiplier) { plotsimple(name, &m_dir, which, multiplier); }
  void plotweighted(std::string &name, int which, int whichredchi2, double multiplier) {
    plotweighted(name, &m_dir, which, whichredchi2, multiplier);
  }
#endif

  // I/O of temporary files for collect mode
  void write(FILE *file, int which = 0);
  void read(FILE *file, int which = 0);

  // these are for the FCN to access what it needs to
  std::vector<double *>::const_iterator residuals_begin() const { return m_residuals.begin(); }
  std::vector<double *>::const_iterator residuals_end() const { return m_residuals.end(); }

  void computeHistogramRangeAndBinning(int which, int &nbins, double &a, double &b);
  void histogramChi2GaussianFit(int which, double &fit_mean, double &fit_sigma);
  void selectPeakResiduals_simple(double nsigma, int nvar, int *vars);
  void selectPeakResiduals(double nsigma, int nvar, int *vars);

  //  void fiducialCuts(double xMin = -1000, double xMax = 1000, double yMin = -1000, double yMax = 1000, bool fidcut1=true);  //  "No fiducial cut"
  //  void fiducialCuts(double xMin = -80.0, double xMax = 80.0, double yMin = -80.0, double yMax = 80.0, bool fidcut1=true);  // "old" fiducial cut
  void fiducialCuts(double xMin = -80.0,
                    double xMax = 80.0,
                    double yMin = -80.0,
                    double yMax = 80.0,
                    bool fidcut1 = false);  // "new" fiducial cut

  virtual void correctBField() = 0;
  virtual void correctBField(int idx_momentum, int idx_q);

  std::vector<bool> &selectedResidualsFlags() { return m_residuals_ok; }

  void eraseNotSelectedResiduals();

protected:
  void initialize_table();
  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;

  int m_residualsModel;
  int m_minHits;
  int m_useResiduals;
  bool m_weightAlignment;
  std::vector<bool> m_fixed;
  std::map<int, double> m_parNum2InitValue;
  int m_printLevel, m_strategy;

  std::vector<double *> m_residuals;
  std::vector<bool> m_residuals_ok;

  std::vector<double> m_value;
  std::vector<double> m_error;
  TMatrixDSym m_cov;
  double m_loglikelihood;

  std::map<int, int> m_parNum2parIdx;

  // center and radii of ellipsoid for peak selection
  // NOTE: 10 is a hardcoded maximum of ellipsoid variables!!!
  // but I can't imagine it ever becoming larger
  double m_center[20];
  double m_radii[20];
};

// A ROOT-sponsored hack to get information into the fit function
class MuonResidualsFitterFitInfo : public TObject {
public:
  MuonResidualsFitterFitInfo(MuonResidualsFitter *fitter) : m_fitter(fitter) {}
  MuonResidualsFitter *fitter() { return m_fitter; }

private:
  MuonResidualsFitter *m_fitter;
#ifdef STANDALONE_FITTER
  ClassDef(MuonResidualsFitterFitInfo, 1);
#endif
};
#ifdef STANDALONE_FITTER
ClassImp(MuonResidualsFitterFitInfo);
#endif

// fit functions (these can't be put in the class; "MuonResidualsFitter_" prefix avoids namespace clashes)
double MuonResidualsFitter_integrate_pureGaussian(double low, double high, double center, double sigma);
double MuonResidualsFitter_logPureGaussian(double residual, double center, double sigma);
Double_t MuonResidualsFitter_pureGaussian_TF1(Double_t *xvec, Double_t *par);
double MuonResidualsFitter_logPureGaussian2D(double x, double y, double x0, double y0, double sx, double sy, double r);
double MuonResidualsFitter_compute_log_convolution(
    double toversigma, double gammaoversigma, double max = 1000., double step = 0.001, double power = 4.);
double MuonResidualsFitter_logPowerLawTails(double residual, double center, double sigma, double gamma);
Double_t MuonResidualsFitter_powerLawTails_TF1(Double_t *xvec, Double_t *par);
double MuonResidualsFitter_logROOTVoigt(double residual, double center, double sigma, double gamma);
Double_t MuonResidualsFitter_ROOTVoigt_TF1(Double_t *xvec, Double_t *par);
double MuonResidualsFitter_logGaussPowerTails(double residual, double center, double sigma);
Double_t MuonResidualsFitter_GaussPowerTails_TF1(Double_t *xvec, Double_t *par);
void MuonResidualsPositionFitter_FCN(int &npar, double *gin, double &fval, double *par, int iflag);
void MuonResidualsAngleFitter_FCN(int &npar, double *gin, double &fval, double *par, int iflag);

#endif  // Alignment_MuonAlignmentAlgorithms_MuonResidualsFitter_H