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:

  static void fcn(int&, double *, double&, double *, int);

  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, 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);
    double initNewN(FitResults);
    double initNewMeanX(FitResults);
    double initNewMeanY(FitResults);
    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);

  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.3
Committed:	Fri Nov 27 03:24:27 2009 UTC (15 years, 5 months ago) by dnisson
Content type:	text/plain
Branch:	MAIN
Changes since 1.2:	+23 -9 lines
Log Message:	Fixed errors in HistoFitter
#	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
17	dnisson	1.3	static void fcn(int&, double , double&, double , int);
18	dnisson	1.2
19	dnisson	1.1	enum TroubleStage {
20			T_NULL = 0,
21			T_SIMPLEX,
22			T_MIGRAD,
23			T_MINOS
24			};
25
26			struct VectorHisto {
27	dnisson	1.2	std::vector< std::vector<double> > bins;
28	dnisson	1.1	double Xlo, Xhi, Ylo, Yhi;
29	dnisson	1.3	};
30	dnisson	1.1
31			struct Trouble {
32			TroubleStage occ;
33			int istat;
34			};
35
36			struct FitResults {
37			TH2 *rebinnedHisto;
38			std::vector< double > chisquare;
39			std::vector< std::vector<TString> > pars;
40			std::vector< std::vector<double> > pval;
41			std::vector< std::vector<double> > perr;
42			std::vector< std::vector<double> > plo;
43			std::vector< std::vector<double> > phi;
44			};
45
46			class ModelDefinition {
47			public:
48			int numberOfGaussians();
49			void setNumberOfGaussians(int);
50			int maxNumberOfGaussians() {
51			return theMaxNumberOfGaussians;
52			}
53			void setMaxNumberOfGaussians(int newNumberOfGaussians) {
54			theMaxNumberOfGaussians = newNumberOfGaussians;
55			}
56
57			double fitFunction(double, double, double *);
58
59	dnisson	1.3	TFormula *indivFormula();
60	dnisson	1.1	void setFormula(TFormula *);
61			void indivParameter(int, std::string&, double&, double&, double&, double&);
62			void setIndivParameter(int, std::string, double, double, double, double);
63	dnisson	1.2	unsigned indivEnergy();
64			void setIndivEnergy(unsigned);
65			unsigned indivMeanX();
66			void setIndivMeanX(unsigned);
67			unsigned indivMeanY();
68			void setIndivMeanY(unsigned);
69	dnisson	1.1	virtual double chisquareSigma(double) = 0;
70			double formulaIntegral(double, double, double, double,
71			double, double, double, double);
72	dnisson	1.3	std::vector< std::vector< double > > subtractCurrentFit(VectorHisto,
73			FitResults);
74	dnisson	1.1	void initParsFromModelDef(std::vector<TString> &,
75			double , double , double , double ,
76			int, FitResults);
77			void findInitialValues(std::vector<TString> &,
78			double , double , double , double ,
79	dnisson	1.3	int, FitResults);
80			double initNewN(FitResults);
81			double initNewMeanX(FitResults);
82			double initNewMeanY(FitResults);
83	dnisson	1.1	void initParameters(TH2 , TMinuit , std::vector<TString> &,
84			double , double , double , double ,
85			int, FitResults);
86
87
88			private:
89			TFormula *theIndivFormula; // formula for individual Gaussian
90			std::vector<std::string> indivParNames;
91			std::vector<double> indivParStarts;
92	dnisson	1.3	std::vector<double> indivParStartingStepSizes;
93	dnisson	1.1	std::vector<double> indivParLowerLimit;
94			std::vector<double> indivParUpperLimit;
95	dnisson	1.2	unsigned theIndivEnergy;
96			unsigned theIndivMeanX;
97			unsigned theIndivMeanY;
98	dnisson	1.1	int theNumberOfGaussians;
99			int theMaxNumberOfGaussians;
100			};
101
102	dnisson	1.2	int getNumberOfDegreesOfFreedom(VectorHisto, int);
103
104	dnisson	1.3	FitResults putParsInFitResults(std::vector<TString>,
105			double , double , double , double );
106			static VectorHisto loadTH2IntoVectorHisto(TH2*);
107			Trouble minimize(TMinuit*, VectorHisto, int);
108
109	dnisson	1.2	FitResults fit_histo(TH2 *hist,
110			std::vector<Trouble> &t,
111			int start_theNumberOfGaussians = 1,
112	dnisson	1.1	double P_cutoff_val = 0.5);
113
114			ModelDefinition& curr_ModelDefinition();
115			void set_ModelDefinition(ModelDefinition *_mdef);
116
117	dnisson	1.2	double fitFunction_TF2(double , double );
118			protected:
119			ModelDefinition *mdef;
120			static ModelDefinition *static_mdef;
121	dnisson	1.3	static VectorHisto energy;
122	dnisson	1.2	};
123	dnisson	1.1
124			#endif