Reference Manual

Go to the documentation of this file.
00001 #ifndef Alignment_MuonAlignmentAlgorithms_MuonResidualsFitter_H
00002 #define Alignment_MuonAlignmentAlgorithms_MuonResidualsFitter_H
00003 
00012 #ifndef STANDALONE_FITTER
00013 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00014 #include "CommonTools/UtilAlgos/interface/TFileService.h"
00015 #include "Alignment/CommonAlignment/interface/Alignable.h"
00016 #endif
00017 
00018 #include "TMinuit.h"
00019 #include "TH1F.h"
00020 #include "TProfile.h"
00021 #include "TF1.h"
00022 #include "TMath.h"
00023 #include "TRandom3.h"
00024 #include "TMatrixDSym.h"
00025 
00026 #include <cstdio>
00027 #include <iostream>
00028 #include <string>
00029 #include <sstream>
00030 #include <map>
00031 
00032 // some mock classes for the case if we want to compile fitters with pure root outside of CMSSW
00033 #ifdef STANDALONE_FITTER
00034 
00035 #include <cassert>
00036 
00037 class Alignable
00038 {
00039 public:
00040   struct Surface
00041   {
00042     double width() {return 300.;}
00043     double length() {return 300.;}
00044   };
00045   Surface s;
00046   Surface & surface() {return s;}
00047 };
00048 
00049 class TFileDirectory
00050 {
00051 public:
00052   template<typename T, typename A1, typename A2, typename A3, typename A4, typename A5>
00053   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; }
00054   template<typename T, typename A1, typename A2, typename A3, typename A4, typename A5, typename A6, typename A7>
00055   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; }
00056   template<typename T, typename A1, typename A2, typename A3, typename A4, typename A5, typename A6, typename A7, typename A8>
00057   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; }
00058 };
00059 
00060 #include <exception>
00061 namespace cms{
00062 class Exception : public std::exception
00063 {
00064 public:
00065   Exception(const std::string & s): name(s) {}
00066   ~Exception () throw () {}
00067   std::string name;
00068   template <class T> Exception& operator<<(const T& what)
00069   {
00070     std::cout<<name<<" exception: "<<what<<std::endl;
00071     return *this;
00072   }
00073 };
00074 }// namespace cms
00075 
00076 #endif // STANDALONE_FITTER
00077 
00078 
00079 
00080 class MuonResidualsFitter
00081 {
00082 public:
00083   enum {
00084     kPureGaussian,
00085     kPowerLawTails,
00086     kROOTVoigt,
00087     kGaussPowerTails,
00088     kPureGaussian2D
00089   };
00090 
00091   enum {
00092     k1DOF,
00093     k5DOF,
00094     k6DOF,
00095     k6DOFrphi,
00096     kPositionFitter,
00097     kAngleFitter,
00098     kAngleBfieldFitter
00099   };
00100 
00101   enum {
00102     k1111,
00103     k1110,
00104     k1100,
00105     k1010,
00106     k0010
00107   };
00108 
00109   struct MuonAlignmentTreeRow
00110   {
00111     Bool_t is_plus; // for +- endcap
00112     Bool_t is_dt; // DT or CSC
00113     UChar_t station; // 8bit uint
00114     Char_t ring_wheel;
00115     UChar_t sector;
00116     Float_t res_x;
00117     Float_t res_y;
00118     Float_t res_slope_x;
00119     Float_t res_slope_y;
00120     Float_t pos_x;
00121     Float_t pos_y;
00122     Float_t angle_x;
00123     Float_t angle_y;
00124     Float_t pz;
00125     Float_t pt;
00126     Char_t q;
00127     Bool_t select;
00128   };
00129 
00130   MuonResidualsFitter(int residualsModel, int minHits, int useResiduals, bool weightAlignment=true)
00131     : m_residualsModel(residualsModel), m_minHits(minHits), m_useResiduals(useResiduals), m_weightAlignment(weightAlignment), m_printLevel(0), m_strategy(1), m_cov(1), m_loglikelihood(0.)
00132   {
00133     if (m_residualsModel != kPureGaussian  &&  m_residualsModel != kPowerLawTails  &&  
00134         m_residualsModel != kROOTVoigt     &&  m_residualsModel != kGaussPowerTails && m_residualsModel != kPureGaussian2D)
00135       throw cms::Exception("MuonResidualsFitter") << "unrecognized residualsModel";
00136   };
00137 
00138   virtual ~MuonResidualsFitter()
00139   {
00140     for (std::vector<double*>::const_iterator residual = residuals_begin();  residual != residuals_end();  ++residual) {
00141       delete [] (*residual);
00142     }
00143   }
00144 
00145   virtual int type() const = 0;
00146   virtual int npar() = 0;
00147   virtual int ndata() = 0;
00148 
00149   int useRes() const { return m_useResiduals; }
00150   int residualsModel() const { return m_residualsModel; }
00151   long numResiduals() const { return m_residuals.size(); }
00152 
00153   void fix(int parNum, bool val=true)
00154   {
00155     assert(0 <= parNum  &&  parNum < npar());
00156     if (m_fixed.size() == 0) m_fixed.resize(npar(), false);
00157     m_fixed[parNum] = val;
00158   }
00159 
00160   bool fixed(int parNum)
00161   {
00162     assert(0 <= parNum  &&  parNum < npar());
00163     if (m_fixed.size() == 0) return false;
00164     else return m_fixed[parNum];
00165   }
00166   int nfixed() { return std::count(m_fixed.begin(), m_fixed.end(), true); }
00167 
00168   void setPrintLevel(int printLevel) { m_printLevel = printLevel; }
00169   void setStrategy(int strategy) { m_strategy = strategy; }
00170 
00171   // an array of the actual residual and associated baggage (qoverpt, trackangle, trackposition)
00172   // arrays passed to fill() are "owned" by MuonResidualsFitter: MuonResidualsFitter will delete them, don't do it yourself!
00173   void fill(double *residual)
00174   {
00175     m_residuals.push_back(residual);
00176     m_residuals_ok.push_back(true);
00177   }
00178 
00179   // this block of results is only valid if fit() returns true
00180   // also gamma is only valid if the model is kPowerLawTails or kROOTVoigt
00181   virtual bool fit(Alignable *ali) = 0;
00182 
00183   double value(int parNum) { assert(0 <= parNum  &&  parNum < npar());  return m_value[parNum]; }
00184   double errorerror(int parNum) { assert(0 <= parNum  &&  parNum < npar());  return m_error[parNum]; }
00185 
00186   // parNum corresponds to the parameter number defined by enums in specific fitters
00187   // parIdx is a continuous index in a set of parameters
00188   int parNum2parIdx(int parNum) { return m_parNum2parIdx[parNum];}
00189 
00190   TMatrixDSym covarianceMatrix() {return m_cov;}
00191   double covarianceElement(int parNum1, int parNum2)
00192   {
00193     assert(0 <= parNum1  &&  parNum1 < npar());
00194     assert(0 <= parNum2  &&  parNum2 < npar());
00195     assert(m_cov.GetNcols() == npar()); // m_cov might have not yet been resized to account for proper #parameters
00196     return m_cov(parNum2parIdx(parNum1),  parNum2parIdx(parNum2));
00197   }
00198   TMatrixDSym correlationMatrix();
00199 
00200   double loglikelihood() { return m_loglikelihood; }
00201 
00202   long numsegments()
00203   {
00204     long num = 0;
00205     for (std::vector<double*>::const_iterator resiter = residuals_begin();  resiter != residuals_end();  ++resiter) num++;
00206     return num;
00207   }
00208 
00209   virtual double sumofweights() = 0;
00210 
00211   // demonstration plots; return reduced chi**2
00212   virtual double plot(std::string name, TFileDirectory *dir, Alignable *ali) = 0;
00213 
00214   void plotsimple(std::string name, TFileDirectory *dir, int which, double multiplier);
00215   void plotweighted(std::string name, TFileDirectory *dir, int which, int whichredchi2, double multiplier);
00216 
00217 #ifdef STANDALONE_FITTER
00218   Alignable m_ali;
00219   TFileDirectory m_dir;
00220   bool fit() { return fit(&m_ali); }
00221   virtual double plot(std::string &name) { return plot(name, &m_dir, &m_ali); }
00222   void plotsimple(std::string &name, int which, double multiplier) { plotsimple(name, &m_dir, which, multiplier); }
00223   void plotweighted(std::string &name, int which, int whichredchi2, double multiplier) { plotweighted(name, &m_dir, which, whichredchi2, multiplier); }
00224 #endif
00225 
00226   // I/O of temporary files for collect mode
00227   void write(FILE *file, int which=0);
00228   void read(FILE *file, int which=0);
00229 
00230   // these are for the FCN to access what it needs to
00231   std::vector<double*>::const_iterator residuals_begin() const { return m_residuals.begin(); }
00232   std::vector<double*>::const_iterator residuals_end() const { return m_residuals.end(); }
00233 
00234   void computeHistogramRangeAndBinning(int which, int &nbins, double &a, double &b); 
00235   void histogramChi2GaussianFit(int which, double &fit_mean, double &fit_sigma);
00236   void selectPeakResiduals_simple(double nsigma, int nvar, int *vars);
00237   void selectPeakResiduals(double nsigma, int nvar, int *vars);
00238 
00239   virtual void correctBField() = 0;
00240   virtual void correctBField(int idx_momentum, int idx_q);
00241 
00242   std::vector<bool> & selectedResidualsFlags() {return m_residuals_ok;}
00243 
00244   void eraseNotSelectedResiduals();
00245 
00246 protected:
00247   void initialize_table();
00248   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);
00249   virtual void inform(TMinuit *tMinuit) = 0;
00250 
00251   int m_residualsModel;
00252   int m_minHits;
00253   int m_useResiduals;
00254   bool m_weightAlignment;
00255   std::vector<bool> m_fixed;
00256   int m_printLevel, m_strategy;
00257 
00258   std::vector<double*> m_residuals;
00259   std::vector<bool> m_residuals_ok;
00260 
00261   std::vector<double> m_value;
00262   std::vector<double> m_error;
00263   TMatrixDSym m_cov;
00264   double m_loglikelihood;
00265 
00266   std::map<int,int> m_parNum2parIdx;
00267 
00268   // center and radii of ellipsoid for peak selection
00269   // NOTE: 10 is a hardcoded maximum of ellipsoid variables!!!
00270   // but I can't imagine it ever becoming larger
00271   double m_center[20];
00272   double m_radii[20];
00273 };
00274 
00275 // A ROOT-sponsored hack to get information into the fit function
00276 class MuonResidualsFitterFitInfo: public TObject
00277 {
00278 public:
00279   MuonResidualsFitterFitInfo(MuonResidualsFitter *fitter) : m_fitter(fitter) {}
00280   MuonResidualsFitter *fitter() { return m_fitter; }
00281 private:
00282   MuonResidualsFitter *m_fitter;
00283 #ifdef STANDALONE_FITTER
00284   ClassDef(MuonResidualsFitterFitInfo,1);
00285 #endif
00286 };
00287 #ifdef STANDALONE_FITTER
00288 ClassImp(MuonResidualsFitterFitInfo);
00289 #endif
00290 
00291 // fit functions (these can't be put in the class; "MuonResidualsFitter_" prefix avoids namespace clashes)
00292 double MuonResidualsFitter_integrate_pureGaussian(double low, double high, double center, double sigma);
00293 double MuonResidualsFitter_logPureGaussian(double residual, double center, double sigma);
00294 Double_t MuonResidualsFitter_pureGaussian_TF1(Double_t *xvec, Double_t *par);
00295 double MuonResidualsFitter_logPureGaussian2D(double x, double y, double x0, double y0, double sx, double sy, double r);
00296 double MuonResidualsFitter_compute_log_convolution(double toversigma, double gammaoversigma, double max=1000., double step=0.001, double power=4.);
00297 double MuonResidualsFitter_logPowerLawTails(double residual, double center, double sigma, double gamma);
00298 Double_t MuonResidualsFitter_powerLawTails_TF1(Double_t *xvec, Double_t *par);
00299 double MuonResidualsFitter_logROOTVoigt(double residual, double center, double sigma, double gamma);
00300 Double_t MuonResidualsFitter_ROOTVoigt_TF1(Double_t *xvec, Double_t *par);
00301 double MuonResidualsFitter_logGaussPowerTails(double residual, double center, double sigma);
00302 Double_t MuonResidualsFitter_GaussPowerTails_TF1(Double_t *xvec, Double_t *par);
00303 void MuonResidualsPositionFitter_FCN(int &npar, double *gin, double &fval, double *par, int iflag);
00304 void MuonResidualsAngleFitter_FCN(int &npar, double *gin, double &fval, double *par, int iflag);
00305 
00306 #endif // Alignment_MuonAlignmentAlgorithms_MuonResidualsFitter_H
/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_5_3_14/src/Alignment/MuonAlignmentAlgorithms/interface/MuonResidualsFitter.h
