| 1 |
|
// Original Author: David Nisson |
| 2 |
< |
// A histogram fitter |
| 2 |
> |
// Fits a two-dimensional model in eta and phi consisting of summed |
| 3 |
> |
// Gaussians to a histogram, supplied in the function fit_histo. |
| 4 |
> |
// This histogram can be of mass, energy, or momentum. |
| 5 |
|
|
| 6 |
|
#ifndef UserCode_JetFitAnalyzer_HistoFitter_h |
| 7 |
|
#define UserCode_JetFitAnalyzer_HistoFitter_h |
| 8 |
|
|
| 9 |
+ |
#include <TH2.h> |
| 10 |
+ |
#include <vector> |
| 11 |
+ |
#include <string> |
| 12 |
+ |
#include <TString.h> |
| 13 |
+ |
#include <TFormula.h> |
| 14 |
+ |
#include <TMinuit.h> |
| 15 |
+ |
|
| 16 |
|
class HistoFitter { |
| 17 |
|
public: |
| 18 |
+ |
// the objective function that gets minimized by TMinuit |
| 19 |
+ |
static void fcn(int&, double *, double&, double *, int); |
| 20 |
+ |
// the correction to Gaussians drifting outside of histogram |
| 21 |
+ |
static double outsiderCorrection(double *, double XbinSize, |
| 22 |
+ |
double YbinSize, |
| 23 |
+ |
double cutoffRad = 0.2); |
| 24 |
|
|
| 25 |
+ |
// place where fitting trouble happened |
| 26 |
|
enum TroubleStage { |
| 27 |
|
T_NULL = 0, |
| 28 |
|
T_SIMPLEX, |
| 30 |
|
T_MINOS |
| 31 |
|
}; |
| 32 |
|
|
| 33 |
+ |
// the histogram class used internally by HistoFitter |
| 34 |
|
struct VectorHisto { |
| 35 |
< |
vector< vector<double> > bins; |
| 35 |
> |
std::vector< std::vector<double> > bins; |
| 36 |
|
double Xlo, Xhi, Ylo, Yhi; |
| 37 |
< |
} energy; |
| 37 |
> |
}; |
| 38 |
|
|
| 39 |
+ |
// information about fitting trouble |
| 40 |
|
struct Trouble { |
| 41 |
|
TroubleStage occ; |
| 42 |
|
int istat; |
| 43 |
|
}; |
| 44 |
|
|
| 45 |
+ |
// structure that stores the fit results |
| 46 |
|
struct FitResults { |
| 47 |
|
TH2 *rebinnedHisto; |
| 48 |
|
std::vector< double > chisquare; |
| 66 |
|
|
| 67 |
|
double fitFunction(double, double, double *); |
| 68 |
|
|
| 69 |
< |
TFormula *getFormula(); |
| 69 |
> |
TFormula *indivFormula(); |
| 70 |
|
void setFormula(TFormula *); |
| 71 |
|
void indivParameter(int, std::string&, double&, double&, double&, double&); |
| 72 |
|
void setIndivParameter(int, std::string, double, double, double, double); |
| 73 |
< |
int indivEnergy(); |
| 74 |
< |
void setIndivEnergy(int); |
| 75 |
< |
int indivMeanX(); |
| 76 |
< |
void setIndivMeanX(int); |
| 77 |
< |
int indivMeanY(); |
| 78 |
< |
void setIndivMeanY(int); |
| 73 |
> |
unsigned indivEnergy(); |
| 74 |
> |
void setIndivEnergy(unsigned); |
| 75 |
> |
unsigned indivMeanX(); |
| 76 |
> |
void setIndivMeanX(unsigned); |
| 77 |
> |
unsigned indivMeanY(); |
| 78 |
> |
void setIndivMeanY(unsigned); |
| 79 |
|
virtual double chisquareSigma(double) = 0; |
| 80 |
|
double formulaIntegral(double, double, double, double, |
| 81 |
< |
double*, double, double, double*); |
| 82 |
< |
std::vector< std::vector< double > > subtractCurrentFit(VectorHisto); |
| 81 |
> |
double*, double, double); |
| 82 |
> |
std::vector< std::vector< double > > subtractCurrentFit(VectorHisto, |
| 83 |
> |
FitResults); |
| 84 |
|
void initParsFromModelDef(std::vector<TString> &, |
| 85 |
|
double *, double *, double *, double *, |
| 86 |
|
int, FitResults); |
| 87 |
|
void findInitialValues(std::vector<TString> &, |
| 88 |
|
double *, double *, double *, double *, |
| 89 |
< |
int); |
| 90 |
< |
double initNewN(); |
| 91 |
< |
double initNewMeanX(); |
| 92 |
< |
double initNewMeanY(); |
| 89 |
> |
int, FitResults, int); |
| 90 |
> |
double initNewN(FitResults, int); |
| 91 |
> |
double initNewMeanX(FitResults, int); |
| 92 |
> |
double initNewMeanY(FitResults, int); |
| 93 |
|
void initParameters(TH2 *, TMinuit *, std::vector<TString> &, |
| 94 |
|
double *, double *, double *, double *, |
| 95 |
|
int, FitResults); |
| 99 |
|
TFormula *theIndivFormula; // formula for individual Gaussian |
| 100 |
|
std::vector<std::string> indivParNames; |
| 101 |
|
std::vector<double> indivParStarts; |
| 102 |
< |
std::vector<double> indivParStartSteps; |
| 102 |
> |
std::vector<double> indivParStartingStepSizes; |
| 103 |
|
std::vector<double> indivParLowerLimit; |
| 84 |
– |
lo; |
| 104 |
|
std::vector<double> indivParUpperLimit; |
| 105 |
< |
int theIndivEnergy; |
| 106 |
< |
int theIndivMeanX; |
| 107 |
< |
int theIndivMeanY; |
| 105 |
> |
unsigned theIndivEnergy; |
| 106 |
> |
unsigned theIndivMeanX; |
| 107 |
> |
unsigned theIndivMeanY; |
| 108 |
|
int theNumberOfGaussians; |
| 109 |
|
int theMaxNumberOfGaussians; |
| 110 |
|
}; |
| 111 |
|
|
| 112 |
< |
FitResults fit_histo(TH2 *hist, std::vector<Trouble> &t, |
| 113 |
< |
void (*cc_minuit)(TMinuit *, TH2 *, int) = 0, |
| 114 |
< |
int start_theNumberOfGaussians = 0, |
| 115 |
< |
int rebinX = 1, int rebinY = 1, |
| 112 |
> |
int getNumberOfDegreesOfFreedom(VectorHisto, int); |
| 113 |
> |
|
| 114 |
> |
FitResults putParsInFitResults(std::vector<TString>, |
| 115 |
> |
double *, double *, double *, double *); |
| 116 |
> |
static VectorHisto loadTH2IntoVectorHisto(TH2*); |
| 117 |
> |
Trouble minimize(TMinuit*, VectorHisto, int); |
| 118 |
> |
|
| 119 |
> |
int findNnewGauss(TH2 *, VectorHisto, FitResults); |
| 120 |
> |
FitResults fit_histo(TH2 *hist, |
| 121 |
> |
std::vector<Trouble> &t, |
| 122 |
> |
int start_theNumberOfGaussians = 1, |
| 123 |
|
double P_cutoff_val = 0.5); |
| 124 |
|
|
| 125 |
|
ModelDefinition& curr_ModelDefinition(); |
| 126 |
|
void set_ModelDefinition(ModelDefinition *_mdef); |
| 127 |
|
|
| 128 |
< |
bool get_ignorezero(); |
| 129 |
< |
void set_ignorezero(bool); |
| 130 |
< |
|
| 131 |
< |
double fit_fcn_TF2(double *, double *); |
| 132 |
< |
// } |
| 107 |
< |
} |
| 128 |
> |
double fitFunction_TF2(double *, double *); |
| 129 |
> |
protected: |
| 130 |
> |
static ModelDefinition *static_mdef; |
| 131 |
> |
static VectorHisto energy; |
| 132 |
> |
}; |
| 133 |
|
|
| 134 |
|
#endif |
