Crombie Tools
PlotFitParameters.cc
Go to the documentation of this file.
1 #include <iostream>
2 #include "TFitResultPtr.h"
3 #include "TFitResult.h"
4 #include "TMath.h"
5 #include "TH2D.h"
6 #include "TH1D.h"
7 #include "TH1F.h"
8 #include "TMatrixDSym.h"
9 #include "TProfile.h"
10 
11 #include "PlotFitParameters.h"
12 
13 ClassImp(PlotFitParameters)
14 
15 //--------------------------------------------------------------------
16 PlotFitParameters::PlotFitParameters() :
17  fFitXBins(0)
18 {
19  fMeans.resize(0);
20  fFunctionComponents.resize(0);
21 }
22 
23 //--------------------------------------------------------------------
24 PlotFitParameters::~PlotFitParameters()
25 {
26  ClearFits();
27 }
28 
29 //--------------------------------------------------------------------
30 void
31 PlotFitParameters::ClearFits()
32 {
33  DisplayFunc(__func__);
34 
35  fMeans.resize(0);
36  fFits.resize(0);
37  fCovs.resize(0);
38 }
39 
40 //--------------------------------------------------------------------
41 void
42 PlotFitParameters::GetMeans(Int_t NumXBins, const Double_t *XBins)
43 {
44  DisplayFunc(__func__);
45 
46  fFitXBins = NumXBins;
47  UInt_t NumPlots = 0;
48 
49  if (fInTrees.size() > 0)
50  NumPlots = fInTrees.size();
51  else if (fInCuts.size() > 0)
52  NumPlots = fInCuts.size();
53  else
54  NumPlots = fInExpr.size();
55 
56  TTree *inTree = fDefaultTree;
57  TString inExpr = fDefaultExpr;
58  TCut inCut = fDefaultCut;
59 
60  for (UInt_t iPlot = 0; iPlot != NumPlots; ++iPlot) {
61  TString tempName;
62  tempName.Form("Hist_%d", fPlotCounter);
63  fPlotCounter++;
64 
65  if (fInTrees.size() != 0)
66  inTree = fInTrees[iPlot];
67  if (fInCuts.size() != 0)
68  inCut = fInCuts[iPlot];
69  if (fInExprXs.size() != 0)
70  fInExprX = fInExprXs[iPlot];
71 
72  TH1 *tempProfile;
73 
74  if (fCutStyle == kBinned) {
75  tempProfile = new TProfile(tempName + "prof", tempName + "prof", NumXBins, XBins);
76  inTree->Draw(fInExprX + ":" + fInExprX + ">>" + tempName + "prof", inCut);
77  }
78  else {
79  tempProfile = new TH1F(tempName + "prof", tempName + "prof", NumXBins, XBins);
80  Int_t addBin = (fCutStyle == kLessThan) ? 1 : 0;
81  for (Int_t iBin = 0; iBin != NumXBins; ++iBin) {
82  tempProfile->SetBinContent(iBin + 1, XBins[iBin + addBin]);
83  tempProfile->SetBinError(iBin + 1, 0);
84  }
85  }
86 
87  fMeans.push_back(tempProfile);
88  }
89 }
90 
91 //--------------------------------------------------------------------
92 void
93 PlotFitParameters::DoFit(TF1* fitFunc, TF1* looseFunc, TH2D* histToFit,
94  TF1**& fitHolder, TMatrixDSym**& covHolder)
95 {
96  DisplayFunc(__func__);
97 
98  Message(eDebug, "Going to fit hist with %f events", histToFit->GetEntries());
99 
100  Int_t NumXBins = histToFit->GetXaxis()->GetNbins();
101  const Double_t *XBins = histToFit->GetXaxis()->GetXbins()->GetArray();
102 
103  TString tempName = histToFit->GetName();
104 
105  if (fMeans.size() == 0)
106  GetMeans(NumXBins, XBins);
107 
108  Double_t MinY = histToFit->GetYaxis()->GetBinLowEdge(1);
109  Double_t MaxY = histToFit->GetYaxis()->GetBinUpEdge(histToFit->GetYaxis()->GetNbins());
110 
111  fitHolder = new TF1*[NumXBins];
112  covHolder = new TMatrixDSym*[NumXBins];
113 
114  for (Int_t iXBin = 0; iXBin != NumXBins; ++iXBin) {
115  Message(eDebug, "Going to fit bin %i", iXBin);
116 
117  TCanvas *tempCanvas = new TCanvas();
118  for (UInt_t iParam = 0; iParam != fGuessParams.size(); ++iParam)
119  fitFunc->SetParameter(fGuessParams[iParam], fGuesses[iParam]);
120 
121  Int_t lowerBin;
122  Int_t upperBin;
123  switch (fCutStyle)
124  {
125  case kBinned:
126  lowerBin = iXBin + 1;
127  upperBin = iXBin + 1;
128  break;
129  case kLessThan:
130  lowerBin = 1;
131  upperBin = iXBin + 1;
132  break;
133  case kGreaterThan:
134  lowerBin = iXBin + 1;
135  upperBin = NumXBins;
136  break;
137  default:
138  std::cout << "What case is that?" << std::endl;
139  exit(1);
140  }
141  TH1D *projection = histToFit->ProjectionY(tempName+"_py", lowerBin, upperBin);
142 
143  Message(eDebug, "Projection has integral %f", projection->GetEntries());
144 
145  if (fLooseFunction != "") {
146  Message(eDebug, "Using loose function: %s", looseFunc->GetTitle());
147 
148  projection->Fit(looseFunc, fFitOptions, "", MinY, MaxY);
149  MapTo(fitFunc, looseFunc);
150  }
151  TFitResultPtr fitResult = projection->Fit(fitFunc, fFitOptions + "S", "", MinY, MaxY);
152  if (fDumpingFits) {
153  TString dumpName;
154  Double_t lower = XBins[lowerBin - 1];
155  Double_t upper = XBins[upperBin];
156  dumpName.Form("DumpFit_%04d_%.2fTo%.2f", fNumFitDumps, lower, upper);
157  std::vector<TF1*> components;
158  for (UInt_t iFunc = 0; iFunc < fFunctionComponents.size(); iFunc++) {
159  TF1 *tempComponent = new TF1(dumpName,fFunctionComponents[iFunc],MinY,MaxY);
160  for (Int_t iParam = 0; iParam != tempComponent->GetNpar(); ++iParam)
161  tempComponent->SetParameter(tempComponent->GetParName(iParam),
162  fitFunc->GetParameter(tempComponent->GetParName(iParam)));
163  tempComponent->Draw("SAME");
164  components.push_back(tempComponent);
165  }
166 
167  dumpName = AddOutDir(dumpName);
168 
169  if (bC)
170  tempCanvas->SaveAs(dumpName+".C");
171  if (bPNG)
172  tempCanvas->SaveAs(dumpName+".png");
173  if (bPDF)
174  tempCanvas->SaveAs(dumpName+".pdf");
175  for (UInt_t iComp = 0; iComp < components.size(); iComp++)
176  delete components[iComp];
177  components.clear();
178  fNumFitDumps++;
179  }
180  fitHolder[iXBin] = (TF1*) fitFunc->Clone();
181  covHolder[iXBin] = new TMatrixDSym(fitFunc->GetNpar());
182  *covHolder[iXBin] = fitResult->GetCovarianceMatrix();
183  }
184 }
185 
186 //--------------------------------------------------------------------
187 std::vector<TGraphErrors*>
188 PlotFitParameters::MakeGraphs(TString ParameterExpr)
189 {
190  DisplayFunc(__func__);
191 
192  Double_t epsilon = 0.001;
193 
194  TGraphErrors *tempGraph;
195  std::vector<TGraphErrors*> theGraphs;
196 
197  TF1 parameterHolder("parameterHolder",ParameterExpr);
198 
199  for (UInt_t iLine = 0; iLine != fFits.size(); ++iLine) {
200  tempGraph = new TGraphErrors(fFitXBins);
201 
202  for (Int_t iXBin = 0; iXBin != fFitXBins; ++iXBin) {
203 
204  // First set the center value of the line being drawn
205  for (Int_t iParam = 0; iParam != parameterHolder.GetNpar(); ++iParam) {
206  Int_t parNumInFit = fFits[iLine][iXBin]->GetParNumber(parameterHolder.GetParName(iParam));
207  parameterHolder.SetParameter(iParam,fFits[iLine][iXBin]->GetParameter(parNumInFit));
208  parameterHolder.SetParError(iParam,fFits[iLine][iXBin]->GetParError(parNumInFit));
209  }
210 
211  tempGraph->SetPoint(iXBin, fMeans[iLine]->GetBinContent(iXBin + 1), parameterHolder.Eval(0));
212 
213  // Next set the error
214  Double_t error2 = 0;
215 
216  for (Int_t iParam = 0; iParam != parameterHolder.GetNpar(); ++iParam) {
217  // Just take the derivative for the first order
218  ParameterExpr.ReplaceAll(TString("[") + TString(parameterHolder.GetParName(iParam)).Strip(TString::kLeading,'p') + TString("]"), 'x');
219  TF1 errorFunc("errorFunc",ParameterExpr);
220  for (Int_t jParam = 0; jParam != errorFunc.GetNpar(); ++jParam)
221  errorFunc.SetParameter(jParam,parameterHolder.GetParameter(errorFunc.GetParName(jParam)));
222  Double_t iDerivative = errorFunc.Derivative(parameterHolder.GetParameter(iParam));
223  Double_t error = iDerivative * parameterHolder.GetParError(iParam);
224  error2 += error*error;
225  ParameterExpr.ReplaceAll('x', TString("[") + TString(parameterHolder.GetParName(iParam)).Strip(TString::kLeading,'p') + TString("]"));
226 
227  for (Int_t jParam = iParam + 1; jParam != parameterHolder.GetNpar(); ++jParam) {
228  // For second order, take the other derivative
229  ParameterExpr.ReplaceAll(TString("[") + TString(parameterHolder.GetParName(jParam)).Strip(TString::kLeading,'p') + TString("]"), 'x');
230 
231  TF1 errorFunc2("errorFunc2",ParameterExpr);
232  for (Int_t kParam = 0; kParam != errorFunc2.GetNpar(); ++kParam)
233  errorFunc2.SetParameter(kParam,parameterHolder.GetParameter(errorFunc.GetParName(kParam)));
234  Double_t jDerivative = errorFunc2.Derivative(parameterHolder.GetParameter(jParam));
235 
236  // Fetch the value from the covariance matrix
237  error2 += 2 * iDerivative * jDerivative *
238  (*fCovs[iLine][iXBin])(fFits[iLine][iXBin]->GetParNumber(parameterHolder.GetParName(iParam)),
239  fFits[iLine][iXBin]->GetParNumber(parameterHolder.GetParName(jParam)));
240 
241  ParameterExpr.ReplaceAll('x', TString("[") + TString(parameterHolder.GetParName(jParam)).Strip(TString::kLeading,'p') + TString("]"));
242  }
243  }
244 
245  tempGraph->SetPointError(iXBin, fMeans[iLine]->GetBinError(iXBin + 1), TMath::Sqrt(error2));
246  }
247  theGraphs.push_back(tempGraph);
248  }
249 
250  return theGraphs;
251 }
252 
253 //--------------------------------------------------------------------
254 std::vector<TGraphErrors*>
255 PlotFitParameters::MakeGraphs(Int_t ParameterNum)
256 {
257  TString ParameterExpr = "";
258  ParameterExpr.Form("[%d]", ParameterNum);
259  return MakeGraphs(ParameterExpr);
260 }