-
Notifications
You must be signed in to change notification settings - Fork 6
/
fit.hh
115 lines (91 loc) · 3.68 KB
/
fit.hh
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
#ifndef __FIT_HH__
#define __FIT_HH__
#include <map>
#include <TMinuit.h>
#include <TRandom3.h>
#include <BAT/BCModel.h>
class TGraph;
class TF1;
class TH1D;
class RandomPrior;
typedef double (*integral_ptr_t)(double*, double*, double*);
class Fitter : public BCModel
{
public:
// constructor/destructor
Fitter();
Fitter(TH1D* data, integral_ptr_t functionIntegral, int maxpar=25);
virtual ~Fitter();
// overloaded methods from BCModel
double LogAPrioriProbability(const std::vector<double> ¶meters);
double LogLikelihood(const std::vector<double> ¶meters);
// setters
void setData(TH1D* hist) { data_=hist; }
void setFunctionIntegral(integral_ptr_t fcnintegral) { functionIntegral_=fcnintegral; }
void setRandomSeed(unsigned int seed) { rand_->SetSeed(seed); }
// getters
bool callLimitReached() { return callLimitReached_; }
const std::string getFitStatus() { return std::string(minuit_.fCstatu.Data()); }
int getNCalls() { return minuit_.fNfcn; }
// parameter manipulation
int defineParameter(int parno, const char* name, double value, double error, double lo, double hi, int isNuisance);
int setParameter(int parno, double value);
int setParLimits(int parno, double loLimit, double hiLimit);
int fixParameter(int parno) { return minuit_.FixParameter(parno); }
int floatParameter(int parno) { return minuit_.Release(parno); }
double getParameter(int parno) const;
double getParError(int parno) const;
void getParLimits(int parno, double &loLimit, double &hiLimit) const;
int getParameter(int parno, double ¤tValue, double ¤tError) const { return minuit_.GetParameter(parno, currentValue, currentError); }
int getNumPars(void) { return minuit_.GetNumPars(); }
double* getParameters(void);
void setPOIIndex(int parno) { poiIndex_=parno; return; }
// fitting, pulls, and all that
void doFit(void);
void doFit(double* emat, int ndim);
TH1D* calcPull(const char* name);
TH1D* calcDiff(const char* name);
// other minor tweaks to Minuit
void setStrategy(int s) { strategy_=s; }
double getStrategy(void) const { return strategy_; }
void setPrintLevel(int p) { printlevel_=p; }
double getPrintLevel(void) const { return printlevel_; }
// evaulate the NLL with the current set of parameters
double evalNLL(void);
// write the TMinuit object
void write(const char* name) const { minuit_.Write(name); return; }
// make pseudodata
TH1D* makePseudoData(const char* name, double* parameters=0);
TH1D* makePseudoDataFromMC(const char* name);
// calculate the posterior distribution
TGraph* calculatePosterior(int nSamples, bool useMCMC=false);
// return the error on a histogram bin with N events
static double histError(double N);
// calculate the CLs
double calculateUpperBoundWithCLs(int nSamples, double alpha);
private:
static void nll(int &, double*, double&, double*, int);
static Fitter* theFitter_;
TRandom3* rand_;
TMinuit minuit_;
integral_ptr_t functionIntegral_;
TH1D* data_;
int printlevel_;
int strategy_;
std::map<int, int> parameterIsNuisance_;
std::map<int, RandomPrior*> priors_; // nuisance parameter priors
double *parameters_;
int poiIndex_;
int nSamples_;
int nCalls_;
bool callLimitReached_;
double poiBestFit_;
double poiUserError_;
bool parRangeSet_; // only used for nuisance parameters with uniform priors
bool useMCMC_;
void evaluateForPosterior(double lo, double mid, double hi, double nllNormalization, std::map<double, double>& fcnEval_);
double computeLikelihoodWithSystematics(double poiVal, double nllNormalization);
// calculate the CLs
std::pair<int, int> calculateCLs_(double poiVal, std::vector<double>& CLb, std::vector<double>& CLsb);
};
#endif