00001
00002
00003
00004
00005
00006
00007
00008
00009
00010 #include <TGraphErrors.h>
00011 #include <TF1.h>
00012 #include <TString.h>
00013 #include <TPad.h>
00014 #include <TLegend.h>
00015
00016 #include "BAT/BCLog.h"
00017 #include "BAT/BCDataSet.h"
00018 #include "BAT/BCDataPoint.h"
00019 #include "BAT/BCMath.h"
00020
00021 #include "BAT/BCGraphFitter.h"
00022
00023
00024
00025 BCGraphFitter::BCGraphFitter() : BCModel("GraphFitter")
00026 {
00027 fGraph = 0;
00028 fFitFunction = 0;
00029 fErrorBand = 0;
00030 fGraphFitFunction = 0;
00031
00032 this -> MCMCSetNIterationsRun(2000);
00033
00034 this -> SetFillErrorBand();
00035 }
00036
00037
00038
00039 BCGraphFitter::BCGraphFitter(TGraphErrors * graph, TF1 * func) : BCModel("GraphFitter")
00040 {
00041 fGraph = 0;
00042 fFitFunction = 0;
00043 fErrorBand = 0;
00044 fGraphFitFunction = 0;
00045
00046 this -> MCMCSetNIterationsRun(2000);
00047
00048 this -> SetGraph(graph);
00049 this -> SetFitFunction(func);
00050
00051 this -> SetFillErrorBand();
00052 }
00053
00054
00055
00056 int BCGraphFitter::SetGraph(TGraphErrors * graph)
00057 {
00058 if(!graph)
00059 {
00060 BCLog::Out(BCLog::error,BCLog::error,"BCGraphFitter::SetGraph() : TGraphErrors not created.");
00061 return 0;
00062 }
00063
00064 int npoints = graph -> GetN();
00065 if(!npoints)
00066 {
00067 BCLog::Out(BCLog::error,BCLog::error,"BCGraphFitter::SetGraph() : TGraphErrors is empty.");
00068 return 0;
00069 }
00070 else if(npoints==1)
00071 {
00072 BCLog::Out(BCLog::error,BCLog::error,"BCGraphFitter::SetGraph() : TGraphErrors has only one point. Not able to fit.");
00073 return 0;
00074 }
00075
00076 fGraph = graph;
00077
00078 double * x = fGraph -> GetX();
00079 double * y = fGraph -> GetY();
00080 double * ex = fGraph -> GetEX();
00081 double * ey = fGraph -> GetEY();
00082
00083 if(!ey)
00084 {
00085 BCLog::Out(BCLog::error,BCLog::error,"BCGraphFitter::SetGraph() : TGraphErrors has NO errors set on Y. Not able to fit.");
00086 return 0;
00087 }
00088
00089 BCDataSet * ds = new BCDataSet();
00090
00091
00092
00093 double xmin=x[0];
00094 double xmax=x[0];
00095 double ymin=y[0];
00096 double ymax=y[0];
00097 for (int i = 0; i < npoints; ++i)
00098 {
00099
00100 double errx = ex ? ex[i] : 0.;
00101
00102
00103 BCDataPoint * dp = new BCDataPoint(4);
00104 dp -> SetValue(0, x[i]);
00105 dp -> SetValue(1, y[i]);
00106 dp -> SetValue(2, errx);
00107 dp -> SetValue(3, ey[i]);
00108 ds -> AddDataPoint(dp);
00109
00110 if(x[i]-errx < xmin)
00111 xmin = x[i]-errx;
00112 else if(x[i]+errx > xmax)
00113 xmax = x[i]+errx;
00114
00115 if(y[i] - 5.*ey[i] < ymin)
00116 ymin = y[i] - 5.*ey[i];
00117 else if(y[i] + 5.*ey[i] > ymax)
00118 ymax = y[i] + 5.*ey[i];
00119 }
00120
00121 this -> SetDataSet(ds);
00122
00123
00124 this -> SetDataBoundaries(0, xmin, xmax);
00125 this -> SetDataBoundaries(1, ymin, ymax);
00126
00127 this -> SetFitFunctionIndices(0, 1);
00128
00129 return this -> GetNDataPoints();
00130 }
00131
00132
00133
00134 int BCGraphFitter::SetFitFunction(TF1 * func)
00135 {
00136 if(!func)
00137 {
00138 BCLog::Out(BCLog::error,BCLog::error,"BCGraphFitter::SetFitFunction() : TF1 not created.");
00139 return 0;
00140 }
00141
00142
00143 int npar = func -> GetNpar();
00144 if(!npar)
00145 {
00146 BCLog::Out(BCLog::error,BCLog::error,"BCGraphFitter::SetFitFunction() : TF1 has zero parameters. Not able to fit.");
00147 return 0;
00148 }
00149
00150
00151 fFitFunction = func;
00152
00153
00154 this -> SetName(TString::Format("GraphFitter with %s",fFitFunction->GetName()));
00155
00156
00157 fParameterSet -> clear();
00158
00159
00160 for (int i = 0; i < npar; ++i)
00161 {
00162 double xmin;
00163 double xmax;
00164 fFitFunction -> GetParLimits(i, xmin, xmax);
00165
00166 this -> AddParameter(fFitFunction->GetParName(i), xmin, xmax);
00167 }
00168
00169 return this -> GetNParameters();
00170 }
00171
00172
00173
00174 BCGraphFitter::~BCGraphFitter()
00175 {}
00176
00177
00178
00179 double BCGraphFitter::LogAPrioriProbability(std::vector <double> parameters)
00180 {
00181
00182 double logprob = 0.;
00183 for(unsigned int i=0; i < this -> GetNParameters(); i++)
00184 logprob -= log(this -> GetParameter(i) -> GetRangeWidth());
00185
00186 return logprob;
00187 }
00188
00189
00190
00191 double BCGraphFitter::LogLikelihood(std::vector <double> params)
00192 {
00193
00194 double logl = 0.;
00195
00196
00197
00198
00199
00200
00201
00202 fFitFunction -> SetParameters(¶ms[0]);
00203
00204
00205 for (int i = 0; i < this -> GetNDataPoints(); i++)
00206 {
00207 std::vector <double> x = GetDataPoint(i) -> GetValues();
00208
00209
00210
00211 double y = x[1];
00212 double yerr = x[3];
00213 double yexp = this -> FitFunction(x,params);
00214
00215
00216
00217 logl += BCMath::LogGaus(y, yexp, yerr, true);
00218 }
00219
00220 return logl;
00221 }
00222
00223
00224
00225 double BCGraphFitter::FitFunction(std::vector <double> x, std::vector <double> params)
00226 {
00227
00228 fFitFunction -> SetParameters(¶ms[0]);
00229
00230 return fFitFunction -> Eval(x[0]);
00231 }
00232
00233
00234
00235 int BCGraphFitter::Fit(TGraphErrors * graph, TF1 * func)
00236 {
00237
00238 if (!this -> SetGraph(graph))
00239 return 0;
00240
00241
00242 if (!this -> SetFitFunction(func))
00243 return 0;
00244
00245
00246 BCLog::Out(BCLog::detail,BCLog::detail,
00247 Form("Fitting %d data points with function of %d parameters",GetNDataPoints(),GetNParameters()));
00248 if(GetNDataPoints() <= (int)GetNParameters())
00249 {
00250 BCLog::Out(BCLog::warning,BCLog::warning,
00251 Form("Number of parameters (%d) lower than or equal to number of points (%d)."
00252 , GetNParameters(), GetNDataPoints()));
00253 BCLog::Out(BCLog::warning,BCLog::warning,"Fit doesn't have much meaning.");
00254 }
00255
00256
00257 this -> MarginalizeAll();
00258
00259
00260
00261 this -> FindModeMinuit( this -> GetBestFitParameters(), -1);
00262
00263
00264
00265
00266 this -> GetPvalueFromChi2(this -> GetBestFitParameters(), 3);
00267 this -> GetPvalueFromChi2NDoF(this -> GetBestFitParameters(), 3);
00268
00269
00270 this -> PrintShortFitSummary(1);
00271
00272 return 1;
00273 }
00274
00275
00276
00277 void BCGraphFitter::DrawFit(const char * options, bool flaglegend)
00278 {
00279 if (!fGraph)
00280 {
00281 BCLog::Out(BCLog::error, BCLog::error,"BCGraphFitter::DrawFit() : TGraphErrors not defined.");
00282 return;
00283 }
00284
00285 if (!fFitFunction)
00286 {
00287 BCLog::Out(BCLog::error, BCLog::error,"BCGraphFitter::DrawFit() : Fit function not defined.");
00288 return;
00289 }
00290
00291
00292 TString opt = options;
00293 opt.ToLower();
00294
00295
00296 if(!opt.Contains("same"))
00297 fGraph -> Draw("ap");
00298
00299
00300 fErrorBand = this -> GetErrorBandGraph(0.16, 0.84);
00301 fErrorBand -> Draw("f same");
00302
00303
00304 fGraphFitFunction = this -> GetFitFunctionGraph( this->GetBestFitParameters() );
00305 fGraphFitFunction -> SetLineColor(kRed);
00306 fGraphFitFunction -> SetLineWidth(2);
00307 fGraphFitFunction -> Draw("l same");
00308
00309
00310
00311 fGraph -> Draw("p same");
00312
00313
00314 if (flaglegend)
00315 {
00316 TLegend * legend = new TLegend(0.25, 0.75, 0.55, 0.95);
00317 legend -> SetBorderSize(0);
00318 legend -> SetFillColor(kWhite);
00319 legend -> AddEntry(fGraph, "Data", "P");
00320 legend -> AddEntry(fGraphFitFunction, "Best fit", "L");
00321 legend -> AddEntry(fErrorBand, "Error band", "F");
00322 legend -> Draw();
00323 }
00324
00325 gPad -> RedrawAxis();
00326 }
00327
00328
00329
00330