JetFitAnalyzer/interface/HistoFitter.h

// Original Author: David Nisson
// A histogram fitter

#ifndef UserCode_JetFitAnalyzer_HistoFitter_h
#define UserCode_JetFitAnalyzer_HistoFitter_h

#include <TH2.h>
#include <vector>
#include <string>
#include <TString.h>
#include <TFormula.h>
#include <TMinuit.h>

class HistoFitter {
public:
  // the objective function that gets minimized by TMinuit
  static void fcn(int&, double *, double&, double *, int);
  // the correction to Gaussians drifting outside of histogram
  static double outsiderCorrection(double *, double XbinSize,
                                   double YbinSize,
                                   double cutoffRad = 0.2);

  enum TroubleStage {
    T_NULL = 0,
    T_SIMPLEX,
    T_MIGRAD,
    T_MINOS
  };

  struct VectorHisto {
    std::vector< std::vector<double> > bins;
    double Xlo, Xhi, Ylo, Yhi;
  };

  struct Trouble {
    TroubleStage occ;
    int istat;
  };

  struct FitResults {
    TH2 *rebinnedHisto;
    std::vector< double > chisquare;
    std::vector< std::vector<TString> > pars;
    std::vector< std::vector<double> > pval;
    std::vector< std::vector<double> > perr;
    std::vector< std::vector<double> > plo;
    std::vector< std::vector<double> > phi;
  };

  class ModelDefinition {
  public:
    int numberOfGaussians();
    void setNumberOfGaussians(int);
    int maxNumberOfGaussians() {
      return theMaxNumberOfGaussians;
    }
    void setMaxNumberOfGaussians(int newNumberOfGaussians) {
      theMaxNumberOfGaussians = newNumberOfGaussians;
    }

    double fitFunction(double, double, double *);

    TFormula *indivFormula();
    void setFormula(TFormula *);
    void indivParameter(int, std::string&, double&, double&, double&, double&);
    void setIndivParameter(int, std::string, double, double, double, double);
    unsigned indivEnergy();
    void setIndivEnergy(unsigned);
    unsigned indivMeanX();
    void setIndivMeanX(unsigned);
    unsigned indivMeanY();
    void setIndivMeanY(unsigned);
    virtual double chisquareSigma(double) = 0;
    double formulaIntegral(double, double, double, double,
                       double*, double, double);
    std::vector< std::vector< double > > subtractCurrentFit(VectorHisto,
                                                            FitResults);
    void initParsFromModelDef(std::vector<TString> &,
                              double *, double *, double *, double *,
                              int, FitResults);
    void findInitialValues(std::vector<TString> &,
                           double *, double *, double *, double *,
                           int, FitResults, int);
    double initNewN(FitResults, int);
    double initNewMeanX(FitResults, int);
    double initNewMeanY(FitResults, int);
    void initParameters(TH2 *, TMinuit *, std::vector<TString> &,
                  double *, double *, double *, double *,
                  int, FitResults);


  private:
    TFormula *theIndivFormula; // formula for individual Gaussian
    std::vector<std::string> indivParNames;
    std::vector<double> indivParStarts;
    std::vector<double> indivParStartingStepSizes;
    std::vector<double> indivParLowerLimit;
    std::vector<double> indivParUpperLimit;
    unsigned theIndivEnergy;
    unsigned theIndivMeanX;
    unsigned theIndivMeanY;
    int theNumberOfGaussians;
    int theMaxNumberOfGaussians;
  };

  int getNumberOfDegreesOfFreedom(VectorHisto, int);

  FitResults putParsInFitResults(std::vector<TString>,
                                 double *, double *, double *, double *);
  static VectorHisto loadTH2IntoVectorHisto(TH2*);
  Trouble minimize(TMinuit*, VectorHisto, int);

  int findNnewGauss(TH2 *, VectorHisto, FitResults);
  FitResults fit_histo(TH2 *hist, 
                        std::vector<Trouble> &t, 
                    int start_theNumberOfGaussians = 1,
                    double P_cutoff_val = 0.5);

  ModelDefinition& curr_ModelDefinition();
  void set_ModelDefinition(ModelDefinition *_mdef);

  double fitFunction_TF2(double *, double *); 
 protected:
  ModelDefinition *mdef;
  static ModelDefinition *static_mdef;
  static VectorHisto energy;
};

#endif
Revision:	1.6
Committed:	Sun Jan 10 01:36:29 2010 UTC (15 years, 3 months ago) by dnisson
Content type:	text/plain
Branch:	MAIN
CVS Tags:	V02-00-00
Changes since 1.5:	+5 -4 lines
Log Message:	Merging Spectrum branch back to trunk
#	User	Rev	Content
1	dnisson	1.1	// Original Author: David Nisson
2			// A histogram fitter
3
4			#ifndef UserCode_JetFitAnalyzer_HistoFitter_h
5			#define UserCode_JetFitAnalyzer_HistoFitter_h
6
7	dnisson	1.3	#include <TH2.h>
8			#include <vector>
9			#include <string>
10			#include <TString.h>
11			#include <TFormula.h>
12			#include <TMinuit.h>
13
14	dnisson	1.1	class HistoFitter {
15			public:
16	dnisson	1.4	// the objective function that gets minimized by TMinuit
17	dnisson	1.3	static void fcn(int&, double , double&, double , int);
18	dnisson	1.4	// the correction to Gaussians drifting outside of histogram
19			static double outsiderCorrection(double *, double XbinSize,
20			double YbinSize,
21			double cutoffRad = 0.2);
22	dnisson	1.2
23	dnisson	1.1	enum TroubleStage {
24			T_NULL = 0,
25			T_SIMPLEX,
26			T_MIGRAD,
27			T_MINOS
28			};
29
30			struct VectorHisto {
31	dnisson	1.2	std::vector< std::vector<double> > bins;
32	dnisson	1.1	double Xlo, Xhi, Ylo, Yhi;
33	dnisson	1.3	};
34	dnisson	1.1
35			struct Trouble {
36			TroubleStage occ;
37			int istat;
38			};
39
40			struct FitResults {
41			TH2 *rebinnedHisto;
42			std::vector< double > chisquare;
43			std::vector< std::vector<TString> > pars;
44			std::vector< std::vector<double> > pval;
45			std::vector< std::vector<double> > perr;
46			std::vector< std::vector<double> > plo;
47			std::vector< std::vector<double> > phi;
48			};
49
50			class ModelDefinition {
51			public:
52			int numberOfGaussians();
53			void setNumberOfGaussians(int);
54			int maxNumberOfGaussians() {
55			return theMaxNumberOfGaussians;
56			}
57			void setMaxNumberOfGaussians(int newNumberOfGaussians) {
58			theMaxNumberOfGaussians = newNumberOfGaussians;
59			}
60
61			double fitFunction(double, double, double *);
62
63	dnisson	1.3	TFormula *indivFormula();
64	dnisson	1.1	void setFormula(TFormula *);
65			void indivParameter(int, std::string&, double&, double&, double&, double&);
66			void setIndivParameter(int, std::string, double, double, double, double);
67	dnisson	1.2	unsigned indivEnergy();
68			void setIndivEnergy(unsigned);
69			unsigned indivMeanX();
70			void setIndivMeanX(unsigned);
71			unsigned indivMeanY();
72			void setIndivMeanY(unsigned);
73	dnisson	1.1	virtual double chisquareSigma(double) = 0;
74			double formulaIntegral(double, double, double, double,
75	dnisson	1.4	double*, double, double);
76	dnisson	1.3	std::vector< std::vector< double > > subtractCurrentFit(VectorHisto,
77			FitResults);
78	dnisson	1.1	void initParsFromModelDef(std::vector<TString> &,
79			double , double , double , double ,
80			int, FitResults);
81			void findInitialValues(std::vector<TString> &,
82			double , double , double , double ,
83	dnisson	1.6	int, FitResults, int);
84			double initNewN(FitResults, int);
85			double initNewMeanX(FitResults, int);
86			double initNewMeanY(FitResults, int);
87	dnisson	1.1	void initParameters(TH2 , TMinuit , std::vector<TString> &,
88			double , double , double , double ,
89			int, FitResults);
90
91
92			private:
93			TFormula *theIndivFormula; // formula for individual Gaussian
94			std::vector<std::string> indivParNames;
95			std::vector<double> indivParStarts;
96	dnisson	1.3	std::vector<double> indivParStartingStepSizes;
97	dnisson	1.1	std::vector<double> indivParLowerLimit;
98			std::vector<double> indivParUpperLimit;
99	dnisson	1.2	unsigned theIndivEnergy;
100			unsigned theIndivMeanX;
101			unsigned theIndivMeanY;
102	dnisson	1.1	int theNumberOfGaussians;
103			int theMaxNumberOfGaussians;
104			};
105
106	dnisson	1.2	int getNumberOfDegreesOfFreedom(VectorHisto, int);
107
108	dnisson	1.3	FitResults putParsInFitResults(std::vector<TString>,
109			double , double , double , double );
110			static VectorHisto loadTH2IntoVectorHisto(TH2*);
111			Trouble minimize(TMinuit*, VectorHisto, int);
112
113	dnisson	1.6	int findNnewGauss(TH2 *, VectorHisto, FitResults);
114	dnisson	1.2	FitResults fit_histo(TH2 *hist,
115			std::vector<Trouble> &t,
116			int start_theNumberOfGaussians = 1,
117	dnisson	1.1	double P_cutoff_val = 0.5);
118
119			ModelDefinition& curr_ModelDefinition();
120			void set_ModelDefinition(ModelDefinition *_mdef);
121
122	dnisson	1.2	double fitFunction_TF2(double , double );
123			protected:
124			ModelDefinition *mdef;
125			static ModelDefinition *static_mdef;
126	dnisson	1.3	static VectorHisto energy;
127	dnisson	1.2	};
128	dnisson	1.1
129			#endif