MuonAlignmentAlgorithms/interface/MuonResidualsFitter.h

0001 #ifndef Alignment_MuonAlignmentAlgorithms_MuonResidualsFitter_H
0002 #define Alignment_MuonAlignmentAlgorithms_MuonResidualsFitter_H
0003
0004 /** \class MuonResidualsFitter
0005  *  $Date: 2011/04/15 21:51:13 $
0006  *  $Revision: 1.16 $
0007  *  \author J. Pivarski - Texas A&M University <pivarski@physics.tamu.edu>
0008  *
0009  *  $Id: MuonResidualsFitter.h,v 1.16 2011/04/15 21:51:13 khotilov Exp $
0010  */
0011
0012 #ifndef STANDALONE_FITTER
0013 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0014 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0015 #include "Alignment/CommonAlignment/interface/Alignable.h"
0016 #endif
0017
0018 #include "TMinuit.h"
0019 #include "TH1F.h"
0020 #include "TProfile.h"
0021 #include "TF1.h"
0022 #include "TMath.h"
0023 #include "TRandom3.h"
0024 #include "TMatrixDSym.h"
0025
0026 #include <cstdio>
0027 #include <iostream>
0028 #include <string>
0029 #include <sstream>
0030 #include <map>
0031
0032 // some mock classes for the case if we want to compile fitters with pure root outside of CMSSW
0033 #ifdef STANDALONE_FITTER
0034
0035 #include <cassert>
0036
0037 class Alignable {
0038 public:
0039   struct Surface {
0040     double width() { return 300.; }
0041     double length() { return 300.; }
0042   };
0043   Surface s;
0044   Surface &surface() { return s; }
0045 };
0046
0047 class TFileDirectory {
0048 public:
0049   template <typename T, typename A1, typename A2, typename A3, typename A4, typename A5>
0050   T *make(const A1 &a1, const A2 &a2, const A3 &a3, const A4 &a4, const A5 &a5) const {
0051     T *t = new T(a1, a2, a3, a4, a5);
0052     return t;
0053   }
0054   template <typename T, typename A1, typename A2, typename A3, typename A4, typename A5, typename A6, typename A7>
0055   T *make(const A1 &a1, const A2 &a2, const A3 &a3, const A4 &a4, const A5 &a5, const A6 &a6, const A7 &a7) const {
0056     T *t = new T(a1, a2, a3, a4, a5, a6, a7);
0057     return t;
0058   }
0059   template <typename T, typename A1, typename A2, typename A3, typename A4, typename A5, typename A6, typename A7, typename A8>
0060   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)
0061       const {
0062     T *t = new T(a1, a2, a3, a4, a5, a6, a7, a8);
0063     return t;
0064   }
0065 };
0066
0067 #include <exception>
0068 namespace cms {
0069   class Exception : public std::exception {
0070   public:
0071     Exception(const std::string &s) : name(s) {}
0072     ~Exception() throw() {}
0073     std::string name;
0074     template <class T>
0075     Exception &operator<<(const T &what) {
0076       std::cout << name << " exception: " << what << std::endl;
0077       return *this;
0078     }
0079   };
0080 }  // namespace cms
0081
0082 #endif  // STANDALONE_FITTER
0083
0084 class MuonResidualsFitter {
0085 public:
0086   enum { kPureGaussian, kPowerLawTails, kROOTVoigt, kGaussPowerTails, kPureGaussian2D };
0087
0088   enum { k1DOF, k5DOF, k6DOF, k6DOFrphi, kPositionFitter, kAngleFitter, kAngleBfieldFitter };
0089
0090   enum { k1111, k1110, k1100, k1010, k0010, k1000, k0100 };
0091
0092   struct MuonAlignmentTreeRow {
0093     Bool_t is_plus;   // for +- endcap
0094     Bool_t is_dt;     // DT or CSC
0095     UChar_t station;  // 8bit uint
0096     Char_t ring_wheel;
0097     UChar_t sector;
0098     Float_t res_x;
0099     Float_t res_y;
0100     Float_t res_slope_x;
0101     Float_t res_slope_y;
0102     Float_t pos_x;
0103     Float_t pos_y;
0104     Float_t angle_x;
0105     Float_t angle_y;
0106     Float_t pz;
0107     Float_t pt;
0108     Char_t q;
0109     Bool_t select;
0110   };
0111
0112   MuonResidualsFitter(int residualsModel, int minHits, int useResiduals, bool weightAlignment = true);
0113   virtual ~MuonResidualsFitter();
0114
0115   virtual int type() const = 0;
0116   virtual int npar() = 0;
0117   virtual int ndata() = 0;
0118
0119   int useRes(int pattern = -1) {
0120     if (pattern >= 0)
0121       m_useResiduals = pattern;
0122     return m_useResiduals;
0123   }
0124   int residualsModel() const { return m_residualsModel; }
0125   long numResiduals() const { return m_residuals.size(); }
0126
0127   void fix(int parNum, bool dofix = true);
0128   bool fixed(int parNum);
0129
0130   int nfixed() { return std::count(m_fixed.begin(), m_fixed.end(), true); }
0131
0132   void setPrintLevel(int printLevel) { m_printLevel = printLevel; }
0133   void setStrategy(int strategy) { m_strategy = strategy; }
0134
0135   void setInitialValue(int parNum, double value) { m_parNum2InitValue[parNum] = value; }
0136
0137   // an array of the actual residual and associated baggage (qoverpt, trackangle, trackposition)
0138   // arrays passed to fill() are "owned" by MuonResidualsFitter: MuonResidualsFitter will delete them, don't do it yourself!
0139   void fill(double *residual);
0140
0141   // this block of results is only valid if fit() returns true
0142   // also gamma is only valid if the model is kPowerLawTails or kROOTVoigt
0143   virtual bool fit(Alignable *ali) = 0;
0144
0145   double value(int parNum) {
0146     assert(0 <= parNum && parNum < npar());
0147     return m_value[parNum];
0148   }
0149   double errorerror(int parNum) {
0150     assert(0 <= parNum && parNum < npar());
0151     return m_error[parNum];
0152   }
0153
0154   // parNum corresponds to the parameter number defined by enums in specific fitters
0155   // parIdx is a continuous index in a set of parameters
0156   int parNum2parIdx(int parNum) { return m_parNum2parIdx[parNum]; }
0157
0158   TMatrixDSym covarianceMatrix() { return m_cov; }
0159   double covarianceElement(int parNum1, int parNum2);
0160   TMatrixDSym correlationMatrix();
0161
0162   double loglikelihood() { return m_loglikelihood; }
0163
0164   long numsegments();
0165
0166   virtual double sumofweights() = 0;
0167
0168   // demonstration plots; return reduced chi**2
0169   virtual double plot(std::string name, TFileDirectory *dir, Alignable *ali) = 0;
0170
0171   void plotsimple(std::string name, TFileDirectory *dir, int which, double multiplier);
0172   void plotweighted(std::string name, TFileDirectory *dir, int which, int whichredchi2, double multiplier);
0173
0174 #ifdef STANDALONE_FITTER
0175   Alignable m_ali;
0176   TFileDirectory m_dir;
0177   bool fit() { return fit(&m_ali); }
0178   virtual double plot(std::string &name) { return plot(name, &m_dir, &m_ali); }
0179   void plotsimple(std::string &name, int which, double multiplier) { plotsimple(name, &m_dir, which, multiplier); }
0180   void plotweighted(std::string &name, int which, int whichredchi2, double multiplier) {
0181     plotweighted(name, &m_dir, which, whichredchi2, multiplier);
0182   }
0183 #endif
0184
0185   // I/O of temporary files for collect mode
0186   void write(FILE *file, int which = 0);
0187   void read(FILE *file, int which = 0);
0188
0189   // these are for the FCN to access what it needs to
0190   std::vector<double *>::const_iterator residuals_begin() const { return m_residuals.begin(); }
0191   std::vector<double *>::const_iterator residuals_end() const { return m_residuals.end(); }
0192
0193   void computeHistogramRangeAndBinning(int which, int &nbins, double &a, double &b);
0194   void histogramChi2GaussianFit(int which, double &fit_mean, double &fit_sigma);
0195   void selectPeakResiduals_simple(double nsigma, int nvar, int *vars);
0196   void selectPeakResiduals(double nsigma, int nvar, int *vars);
0197
0198   //  void fiducialCuts(double xMin = -1000, double xMax = 1000, double yMin = -1000, double yMax = 1000, bool fidcut1=true);  //  "No fiducial cut"
0199   //  void fiducialCuts(double xMin = -80.0, double xMax = 80.0, double yMin = -80.0, double yMax = 80.0, bool fidcut1=true);  // "old" fiducial cut
0200   void fiducialCuts(double xMin = -80.0,
0201                     double xMax = 80.0,
0202                     double yMin = -80.0,
0203                     double yMax = 80.0,
0204                     bool fidcut1 = false);  // "new" fiducial cut
0205
0206   virtual void correctBField() = 0;
0207   virtual void correctBField(int idx_momentum, int idx_q);
0208
0209   std::vector<bool> &selectedResidualsFlags() { return m_residuals_ok; }
0210
0211   void eraseNotSelectedResiduals();
0212
0213 protected:
0214   void initialize_table();
0215   bool dofit(void (*fcn)(int &, double *, double &, double *, int),
0216              std::vector<int> &parNum,
0217              std::vector<std::string> &parName,
0218              std::vector<double> &start,
0219              std::vector<double> &step,
0220              std::vector<double> &low,
0221              std::vector<double> &high);
0222   virtual void inform(TMinuit *tMinuit) = 0;
0223
0224   int m_residualsModel;
0225   int m_minHits;
0226   int m_useResiduals;
0227   bool m_weightAlignment;
0228   std::vector<bool> m_fixed;
0229   std::map<int, double> m_parNum2InitValue;
0230   int m_printLevel, m_strategy;
0231
0232   std::vector<double *> m_residuals;
0233   std::vector<bool> m_residuals_ok;
0234
0235   std::vector<double> m_value;
0236   std::vector<double> m_error;
0237   TMatrixDSym m_cov;
0238   double m_loglikelihood;
0239
0240   std::map<int, int> m_parNum2parIdx;
0241
0242   // center and radii of ellipsoid for peak selection
0243   // NOTE: 10 is a hardcoded maximum of ellipsoid variables!!!
0244   // but I can't imagine it ever becoming larger
0245   double m_center[20];
0246   double m_radii[20];
0247 };
0248
0249 // A ROOT-sponsored hack to get information into the fit function
0250 class MuonResidualsFitterFitInfo : public TObject {
0251 public:
0252   MuonResidualsFitterFitInfo(MuonResidualsFitter *fitter) : m_fitter(fitter) {}
0253   MuonResidualsFitter *fitter() { return m_fitter; }
0254
0255 private:
0256   MuonResidualsFitter *m_fitter;
0257 #ifdef STANDALONE_FITTER
0258   ClassDef(MuonResidualsFitterFitInfo, 1);
0259 #endif
0260 };
0261 #ifdef STANDALONE_FITTER
0262 ClassImp(MuonResidualsFitterFitInfo);
0263 #endif
0264
0265 // fit functions (these can't be put in the class; "MuonResidualsFitter_" prefix avoids namespace clashes)
0266 double MuonResidualsFitter_integrate_pureGaussian(double low, double high, double center, double sigma);
0267 double MuonResidualsFitter_logPureGaussian(double residual, double center, double sigma);
0268 Double_t MuonResidualsFitter_pureGaussian_TF1(Double_t *xvec, Double_t *par);
0269 double MuonResidualsFitter_logPureGaussian2D(double x, double y, double x0, double y0, double sx, double sy, double r);
0270 double MuonResidualsFitter_compute_log_convolution(
0271     double toversigma, double gammaoversigma, double max = 1000., double step = 0.001, double power = 4.);
0272 double MuonResidualsFitter_logPowerLawTails(double residual, double center, double sigma, double gamma);
0273 Double_t MuonResidualsFitter_powerLawTails_TF1(Double_t *xvec, Double_t *par);
0274 double MuonResidualsFitter_logROOTVoigt(double residual, double center, double sigma, double gamma);
0275 Double_t MuonResidualsFitter_ROOTVoigt_TF1(Double_t *xvec, Double_t *par);
0276 double MuonResidualsFitter_logGaussPowerTails(double residual, double center, double sigma);
0277 Double_t MuonResidualsFitter_GaussPowerTails_TF1(Double_t *xvec, Double_t *par);
0278 void MuonResidualsPositionFitter_FCN(int &npar, double *gin, double &fval, double *par, int iflag);
0279 void MuonResidualsAngleFitter_FCN(int &npar, double *gin, double &fval, double *par, int iflag);
0280
0281 #endif  // Alignment_MuonAlignmentAlgorithms_MuonResidualsFitter_H