Crombie Tools
Plot2D.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 
7 #include "Plot2D.h"
8 
9 ClassImp(Plot2D)
10 
11 //--------------------------------------------------------------------
12 Plot2D::Plot2D() :
13  fInExprX(""),
14  fFunctionString(""),
15  fLooseFunction(""),
16  fDumpingFits(false),
17  fNumFitDumps(0),
18  fNumPoints(100),
19  fFitOptions("MLEQ")
20 {
21 
22  fFits.resize(0);
23  fCovs.resize(0);
24  fGuessParams.resize(0);
25  fGuesses.resize(0);
26  fLooseGuessParams.resize(0);
27  fLooseGuesses.resize(0);
28  fParamFrom.resize(0);
29  fParamTo.resize(0);
30  fParams.resize(0);
31  fParamLows.resize(0);
32  fParamHighs.resize(0);
33  fLooseParams.resize(0);
34  fLooseParamLows.resize(0);
35  fLooseParamHighs.resize(0);
36  fInExprXs.resize(0);
37 
38 }
39 
40 //--------------------------------------------------------------------
41 Plot2D::~Plot2D()
42 {
43 
44  ClearFits();
45 
46 }
47 
48 //--------------------------------------------------------------------
49 void
50 Plot2D::ClearFits()
51 {
52 
53  fFits.resize(0);
54  fCovs.resize(0);
55 
56 }
57 
58 //--------------------------------------------------------------------
59 void
60 Plot2D::ResetParameterLimits()
61 {
62 
63  fParams.resize(0);
64  fParamLows.resize(0);
65  fParamHighs.resize(0);
66 
67 }
68 
69 //--------------------------------------------------------------------
70 void
71 Plot2D::SetParameterLimits(Int_t param, Double_t low, Double_t high)
72 {
73 
74  DisplayFunc(__func__);
75  fParams.push_back(param);
76  fParamLows.push_back(low);
77  fParamHighs.push_back(high);
78 
79 }
80 
81 //--------------------------------------------------------------------
82 void
83 Plot2D::SetLooseLimits(Int_t param, Double_t low, Double_t high)
84 {
85  DisplayFunc(__func__);
86  fLooseParams.push_back(param);
87  fLooseParamLows.push_back(low);
88  fLooseParamHighs.push_back(high);
89 }
90 
91 //--------------------------------------------------------------------
92 void
93 Plot2D::AddMapping(Int_t from, Int_t to, Double_t fact)
94 {
95 
96  DisplayFunc(__func__);
97  fParamFrom.push_back(from);
98  fParamTo.push_back(to);
99  fParamFactor.push_back(fact);
100 
101 }
102 
103 //--------------------------------------------------------------------
104 void
105 Plot2D::MapTo(TF1* fitFunc, TF1* looseFunc)
106 {
107 
108  DisplayFunc(__func__);
109 
110  for (UInt_t iParam = 0; iParam != fParamFrom.size(); ++iParam) {
111  double value = looseFunc->GetParameter(fParamFrom[iParam]) * fParamFactor[iParam];
112  Message(eDebug, "Mapping param %i times %f to %i; value %f",
113  fParamFrom[iParam], fParamFactor[iParam], fParamTo[iParam], value);
114  fitFunc->SetParameter(fParamTo[iParam], value);
115  }
116 
117 }
118 
119 //--------------------------------------------------------------------
120 void
121 Plot2D::DoFit(TF1* fitFunc, TF1* looseFunc, TH2D* histToFit,
122  TF1**& fitHolder, TMatrixDSym**& covHolder)
123 {
124 
125  DisplayFunc(__func__);
126  fitHolder = new TF1*[1];
127  covHolder = new TMatrixDSym*[1];
128  TCanvas *tempCanvas = new TCanvas();
129 
130  if (fLooseFunction != "") {
131 
132  histToFit->Fit(looseFunc,fFitOptions);
133  MapTo(fitFunc,looseFunc);
134 
135  }
136 
137  TFitResultPtr fitResult = histToFit->Fit(fitFunc,fFitOptions + "S");
138 
139  fitHolder[0] = (TF2*) fitFunc->Clone();
140  covHolder[0] = new TMatrixDSym(fitFunc->GetNpar());
141  *covHolder[0] = fitResult->GetCovarianceMatrix();
142 
143  if (fDumpingFits) {
144 
145  TString dumpName;
146  dumpName.Form("DumpFit_%04d_2D", fNumFitDumps++);
147 
148  dumpName = AddOutDir(dumpName);
149 
150  if (bC)
151  tempCanvas->SaveAs(dumpName+".C");
152  if (bPNG)
153  tempCanvas->SaveAs(dumpName+".png");
154  if (bPDF)
155  tempCanvas->SaveAs(dumpName+".pdf");
156 
157  }
158 
159 }
160 
161 //--------------------------------------------------------------------
162 void
163 Plot2D::DoFits(Int_t NumXBins, Double_t *XBins,
164  Int_t NumYBins, Double_t MinY, Double_t MaxY)
165 {
166 
167  DisplayFunc(__func__);
168  ClearFits();
169  UInt_t NumPlots = 0;
170 
171  if (fFunctionString == "") {
172 
173  std::cerr << "You haven't set a function!" << std::endl;
174  exit(256);
175 
176  }
177  if (fInExprX == "" && fInExprXs.size() == 0) {
178 
179  std::cerr << "You haven't initialized an x expression yet!" << std::endl;
180  exit(257);
181 
182  }
183 
184  if (fInTrees.size() != 0)
185  NumPlots = fInTrees.size();
186  else if (fInCuts.size() != 0)
187  NumPlots = fInCuts.size();
188  else
189  NumPlots = fInExpr.size();
190 
191  if(NumPlots == 0){
192 
193  std::cerr << "Nothing has been initialized in resolution plot." << std::endl;
194  exit(300);
195 
196  }
197 
198  TTree *inTree = fDefaultTree;
199  TString inExpr = fDefaultExpr;
200  TCut inCut = fDefaultCut;
201 
202  TH2D *tempHist;
203  TProfile *tempProfile;
204 
205  TF1 *looseFunc = NULL;
206 
207  if (fLooseFunction != "") {
208 
209  looseFunc = MakeFunction(fLooseFunction, XBins[0], XBins[NumXBins - 1], MinY, MaxY);
210  looseFunc->SetLineColor(kGreen);
211 
212  }
213 
214  for (UInt_t iGuess = 0; iGuess != fLooseGuessParams.size(); ++iGuess)
215  looseFunc->SetParameter(fLooseGuessParams[iGuess], fLooseGuesses[iGuess]);
216 
217  for (UInt_t iParam = 0; iParam != fLooseParams.size(); ++iParam)
218  looseFunc->SetParLimits(fLooseParams[iParam], fLooseParamLows[iParam], fLooseParamHighs[iParam]);
219 
220  TF1 *fitFunc = MakeFunction(fFunctionString, XBins[0], XBins[NumXBins - 1], MinY, MaxY);
221  fitFunc->SetLineColor(kBlue);
222 
223  for (UInt_t iGuess = 0; iGuess != fGuessParams.size(); ++iGuess)
224  fitFunc->SetParameter(fGuessParams[iGuess], fGuesses[iGuess]);
225 
226  for (UInt_t iParam = 0; iParam != fParams.size(); ++iParam)
227  fitFunc->SetParLimits(fParams[iParam], fParamLows[iParam], fParamHighs[iParam]);
228 
229  TF1 **holdFits;
230  TMatrixDSym **holdCovs;
231 
232  for (UInt_t iPlot = 0; iPlot < NumPlots; iPlot++) {
233 
234  std::cout << NumPlots - iPlot << " more to go." << std::endl;
235 
236  if (fInTrees.size() != 0)
237  inTree = fInTrees[iPlot];
238  if (fInCuts.size() != 0)
239  inCut = fInCuts[iPlot];
240  if (fInExpr.size() != 0)
241  inExpr = fInExpr[iPlot];
242  if (fInExprXs.size() != 0)
243  fInExprX = fInExprXs[iPlot];
244 
245  TString tempName;
246  tempName.Form("Hist_%d", fPlotCounter);
247  fPlotCounter++;
248  tempHist = new TH2D(tempName, tempName, NumXBins, XBins, NumYBins, MinY, MaxY);
249  tempHist->Sumw2();
250 
251  TString drawCommand = inExpr + ":" + fInExprX + ">>" + tempName;
252 
253  Message(eDebug, "About to use tree at %p to draw %s with cut %s",
254  inTree, drawCommand.Data(), inCut.GetTitle());
255 
256  inTree->Draw(drawCommand, inCut);
257  Message(eDebug, "Histogram has %f events", tempHist, tempHist->Integral());
258 
259  TString dumpTitle = fLegendEntries[iPlot] + ";" + inExpr + ";Num Events";
260 
261  // if (fLooseFunction != "")
262  // looseFunc->SetTitle(dumpTitle);
263 
264  fitFunc->SetTitle(dumpTitle);
265  tempHist->SetTitle(fLegendEntries[iPlot] + ";" + fInExprX + ";" + inExpr + ";Num Events");
266 
267  DoFit(fitFunc, looseFunc, tempHist, holdFits, holdCovs);
268 
269  fFits.push_back(holdFits);
270  fCovs.push_back(holdCovs);
271  delete tempHist;
272 
273  }
274 
275 }
276 
277 //--------------------------------------------------------------------
278 void
279 Plot2D::DoFits(Int_t NumXBins, Double_t MinX, Double_t MaxX,
280  Int_t NumYBins, Double_t MinY, Double_t MaxY)
281 {
282 
283  Double_t XBins[NumXBins+1];
284  ConvertToArray(NumXBins, MinX, MaxX, XBins);
285  DoFits(NumXBins, XBins, NumYBins, MinY, MaxY);
286 
287 }
288 
289 //--------------------------------------------------------------------
290 std::vector<TF1*>
291 Plot2D::MakeFuncs(TString ParameterExpr, Double_t MinX, Double_t MaxX)
292 {
293 
294  DisplayFunc(__func__);
295  TF1 *tempFunc;
296  std::vector<TF1*> theFuncs;
297 
298  for (UInt_t iLine = 0; iLine != fFits.size(); ++iLine) {
299 
300  tempFunc = new TF1("parameterHolder", ParameterExpr, MinX, MaxX);
301 
302  for (Int_t iParam = 0; iParam != tempFunc->GetNpar(); ++iParam) {
303 
304  Int_t parNumInFit = fFits[iLine][0]->GetParNumber(tempFunc->GetParName(iParam));
305  tempFunc->SetParameter(iParam, fFits[iLine][0]->GetParameter(parNumInFit));
306  tempFunc->SetParError(iParam, fFits[iLine][0]->GetParError(parNumInFit));
307 
308  }
309 
310  theFuncs.push_back(tempFunc);
311 
312  }
313 
314  return theFuncs;
315 
316 }
317 
318 //--------------------------------------------------------------------
319 std::vector<TGraphErrors*>
320 Plot2D::MakeGraphs(TString ParameterExpr)
321 {
322 
323  DisplayFunc(__func__);
324  Double_t MinX = 0.0;
325  Double_t MaxX = 0.0;
326 
327  fFits[0][0]->GetRange(MinX, MaxX);
328 
329  std::vector<TF1*> theFuncs = MakeFuncs(ParameterExpr, MinX, MaxX);
330  std::vector<TGraphErrors*> theGraphs;
331  TGraphErrors* tempGraph;
332 
333  Double_t width = (MaxX - MinX)/fNumPoints;
334 
335  for (UInt_t iFunc = 0; iFunc != theFuncs.size(); ++iFunc) {
336 
337  std::vector<TF1*> holdFuncs;
338 
339  for (Int_t iParam = 0; iParam != theFuncs[iFunc]->GetNpar(); ++iParam) {
340 
341  TString parToReplace = TString('[') + theFuncs[iFunc]->GetParName(iParam) + TString(']');
342  TF1* holdFunc = new TF1("holding",theFuncs[iFunc]->GetExpFormula().ReplaceAll('x',"[x0]").ReplaceAll(parToReplace,'x'));
343 
344  for (Int_t jParam = 0; jParam != holdFunc->GetNpar(); ++jParam) {
345 
346  TString paramName = holdFunc->GetParName(jParam);
347 
348  if (paramName == "x0")
349  continue;
350 
351  holdFunc->SetParameter(jParam,theFuncs[iFunc]->GetParameter(paramName));
352 
353  }
354 
355  holdFuncs.push_back(holdFunc);
356 
357  }
358 
359  tempGraph = new TGraphErrors(fNumPoints);
360  Double_t xVal = MinX;
361 
362  for (Int_t iPoint = 0; iPoint != fNumPoints + 1; ++iPoint) {
363 
364  tempGraph->SetPoint(iPoint,xVal,theFuncs[iFunc]->Eval(xVal));
365  Double_t error2 = 0;
366 
367  for (Int_t iParam = 0; iParam != theFuncs[iFunc]->GetNpar(); ++iParam) {
368 
369  holdFuncs[iParam]->SetParameter("x0",xVal);
370 
371  for (Int_t jParam = 0; jParam != iParam + 1; ++jParam) {
372 
373  Double_t factor = 2.0;
374 
375  if (iParam == jParam)
376  factor = 1.0;
377 
378  Double_t toAdd = holdFuncs[iParam]->Derivative(theFuncs[iFunc]->GetParameter(iParam)) * ((*(fCovs[iFunc][0]))(iParam, jParam));
379  error2 += toAdd*toAdd;
380 
381  }
382 
383  }
384 
385  tempGraph->SetPointError(iPoint,0,TMath::Sqrt(error2));
386  xVal += width;
387 
388  }
389 
390  theGraphs.push_back(tempGraph);
391 
392  for (UInt_t iParam = 0; iParam != holdFuncs.size(); ++iParam)
393  delete holdFuncs[iParam];
394 
395  }
396 
397  for (UInt_t iFunc = 0; iFunc != theFuncs.size(); ++iFunc)
398  delete theFuncs[iFunc];
399 
400  return theGraphs;
401 
402 }
403 
404 //--------------------------------------------------------------------
405 void
406 Plot2D::MakeCanvas(TString FileBase, std::vector<TGraphErrors*> theGraphs, TString XLabel, TString YLabel,
407  Double_t YMin, Double_t YMax, Bool_t logY)
408 {
409 
410  for (UInt_t iGraph = 0; iGraph != theGraphs.size(); ++iGraph)
411  theGraphs[iGraph]->GetYaxis()->SetRangeUser(YMin, YMax);
412 
413  BaseCanvas(FileBase, theGraphs, XLabel, YLabel, logY);
414 
415 }
416 
417 //--------------------------------------------------------------------
418 void
419 Plot2D::MakeCanvas(TString FileBase, TString ParameterExpr, TString XLabel, TString YLabel,
420  Double_t YMin, Double_t YMax, Bool_t logY)
421 {
422 
423  std::vector<TGraphErrors*> theGraphs = MakeGraphs(ParameterExpr);
424  MakeCanvas(FileBase, theGraphs, XLabel, YLabel, YMin, YMax, logY);
425 
426  for (UInt_t iGraph = 0; iGraph != theGraphs.size(); ++iGraph)
427  delete theGraphs[iGraph];
428 
429 }
430 
431 //--------------------------------------------------------------------
432 void
433 Plot2D::MakeCanvas(TString FileBase, Int_t ParameterNum, TString XLabel, TString YLabel,
434  Double_t YMin, Double_t YMax, Bool_t logY)
435 {
436 
437  TString ParameterExpr = "";
438  ParameterExpr.Form("[%d]", ParameterNum);
439  MakeCanvas(FileBase, ParameterExpr, XLabel, YLabel, YMin, YMax, logY);
440 
441 }