source: trunk/MagicSoft/Mars/macros/estfit.C@ 1656

Last change on this file since 1656 was 1633, checked in by tbretz, 23 years ago
*** empty log message ***
File size: 8.4 KB
Line 
1#include <fstream.h>
2
3#include "MTask.h"
4#include "MParList.h"
5
6#include "MDataChain.h"
7
8class MChisqEval : public MTask
9{
10private:
11 static const TString gsDefName;
12 static const TString gsDefTitle;
13
14 Double_t fChisq; //! Evaluated chi square
15
16 MData *fData0; // Data Member one (monte carlo data or chisq function)
17 MData *fData1; // Data Member two (measured data)
18
19 // --------------------------------------------------------------------------
20 //
21 // Implementation of SavePrimitive. Used to write the call to a constructor
22 // to a macro. In the original root implementation it is used to write
23 // gui elements to a macro-file.
24 //
25 void StreamPrimitive(ofstream &out) const
26 {
27 out << " MChisqEval " << GetUniqueName() << ";";
28 if (fData0)
29 out << " " << GetUniqueName() << ".SetY1(\"" << fData0->GetRule() << "\");" << endl;
30 if (fData1)
31 out << " " << GetUniqueName() << ".SetY1(\"" << fData1->GetRule() << "\");" << endl;
32 }
33
34 enum { kIsOwner = BIT(14) };
35
36public:
37 MChisqEval(const char *name=NULL, const char *title=NULL) : fData0(NULL), fData1(NULL)// : fMatrix(mat), fColumn(col), fEvalE(evale)
38 {
39 fName = name ? name : gsDefName.Data();
40 fTitle = title ? title : gsDefTitle.Data();
41 }
42
43 MChisqEval(MData *y1, const char *name=NULL, const char *title=NULL) : fData0(NULL), fData1(NULL)// : fMatrix(mat), fColumn(col), fEvalE(evale)
44 {
45 fName = name ? name : gsDefName.Data();
46 fTitle = title ? title : gsDefTitle.Data();
47 SetY1(y1);
48 }
49 MChisqEval(MData *y1, MData *y2, const char *name=NULL, const char *title=NULL) : fData0(NULL), fData1(NULL)// : fMatrix(mat), fColumn(col), fEvalE(evale)
50 {
51 fName = name ? name : gsDefName.Data();
52 fTitle = title ? title : gsDefTitle.Data();
53 SetY1(y1);
54 SetY2(y2);
55 }
56 ~MChisqEval()
57 {
58 if (fData0 && (fData0->TestBit(kCanDelete) || TestBit(kIsOwner)))
59 delete fData0;
60
61 if (fData1 && (fData1->TestBit(kCanDelete) || TestBit(kIsOwner)))
62 delete fData1;
63 }
64
65 void SetY1(MData *data)
66 {
67 // Set MC value
68 if (fData0 && (fData0->TestBit(kCanDelete) || TestBit(kIsOwner)))
69 delete fData0;
70 fData0 = data;
71 fData0->SetBit(kCanDelete);
72 AddToBranchList(fData0->GetDataMember());
73 }
74 void SetY2(MData *data)
75 {
76 // Set measured/estimated value
77 if (fData1 && (fData1->TestBit(kCanDelete) || TestBit(kIsOwner)))
78 delete fData1;
79 fData1 = data;
80 fData1->SetBit(kCanDelete);
81 AddToBranchList(fData1->GetDataMember());
82 }
83 void SetY1(const TString data) { SetY1(new MDataChain(data)); }
84 void SetY2(const TString data) { SetY2(new MDataChain(data)); }
85
86 void SetOwner(Bool_t o=kTRUE) { o ? SetBit(kIsOwner) : ResetBit(kIsOwner); }
87
88 Bool_t PreProcess(MParList *plist)
89 {
90 fChisq = 0;
91
92 if (!fData0)
93 return kFALSE;
94
95 if (!fData0->PreProcess(plist))
96 return kFALSE;
97
98 if (fData1)
99 if (!fData1->PreProcess(plist))
100 return kFALSE;
101
102 return kTRUE;
103 }
104
105 Bool_t Process()
106 {
107 const Double_t y1 = fData0->GetValue();
108 const Double_t y2 = fData1 ? fData1->GetValue() : 0;
109
110 const Double_t dy = y2-y1;
111 const Double_t err = fData1 ? y1*y1 : 1;
112
113 fChisq += dy*dy/err;
114 return kTRUE;
115 }
116
117 Bool_t PostProcess()
118 {
119 fChisq /= GetNumExecutions();
120 return kTRUE;
121 }
122
123 Double_t GetChisq() const { return fChisq; }
124
125 ClassDef(MChisqEval, 0)
126};
127
128ClassImp(MChisqEval);
129
130const TString MChisqEval::gsDefName = "MChisqEval";
131const TString MChisqEval::gsDefTitle = "Evaluate a chisq";
132
133// --------------------------------------------------------------------------
134
135#include <TMinuit.h>
136
137#include "MEvtLoop.h"
138#include "MTaskList.h"
139#include "MEnergyEstParam.h"
140
141void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
142{
143 MEvtLoop &evtloop = *(MEvtLoop*)gMinuit->GetObjectFit();
144
145 MTaskList &tlist = *(MTaskList*)evtloop.GetParList()->FindObject("MTaskList"); //GetTaskList();
146
147 MChisqEval &eval = *(MChisqEval*) tlist.FindObject("MChisqEval");
148 MEnergyEstParam &eest = *(MEnergyEstParam*)tlist.FindObject("MEnergyEstParam");
149
150 eest.SetCoeff(TArrayD(eest.GetNumCoeff(), par));
151
152 evtloop.Eventloop();
153
154 f = eval.GetChisq();
155}
156
157// --------------------------------------------------------------------------
158#include <TStopwatch.h>
159
160#include "MHMatrix.h"
161#include "MEnergyEst.h"
162
163#include "MDataMember.h"
164#include "MDataElement.h"
165
166#include "MReadTree.h"
167#include "MMatrixLoop.h"
168
169void estfit(Bool_t evalenergy=kFALSE)
170{
171 //
172 // Fill events into a MHMatrix
173 //
174 MParList parlist;
175 MHMatrix matrix;
176
177 Int_t col = matrix.AddColumn(evalenergy?"MMcEvt.fEnergy":"MMcEvt.fImpact");
178
179 MEnergyEstParam eest;
180 eest.Add("MHillasSrc");
181 eest.InitMapping(&matrix);
182
183 MReadTree read("Events", "gammas.root");
184 read.DisableAutoScheme();
185
186 if (!matrix.Fill(&parlist, &read))
187 return;
188
189 //
190 // Setup the tasklist used to evaluate the needed chisq
191 //
192 MTaskList tasklist;
193 parlist.AddToList(&tasklist);
194
195 MMatrixLoop loop(&matrix);
196
197 MChisqEval eval;
198 eval.SetY1(new MDataElement(&matrix, col));
199 eval.SetY2(new MDataMember(evalenergy ? "MEnergyEst.fEnergy" : "MEnergyEst.fImpact"));
200 eval.SetOwner();
201
202 tasklist.AddToList(&loop);
203 tasklist.AddToList(&eest);
204 tasklist.AddToList(&eval);
205
206 MEvtLoop evtloop;
207 evtloop.SetParList(&parlist);
208
209 //
210 // Be carefull: This is not thread safe
211 //
212 TMinuit minuit(12);
213 minuit.SetFCN(fcn);
214
215 // Ready for: minuit.mnexcm("SET ERR", arglist, 1, ierflg)
216 if (minuit.SetErrorDef(1))
217 {
218 cout << "SetErrorDef failed." << endl;
219 return;
220 }
221
222 //
223 // Set initial values
224 //
225 TArrayD fA(5);
226 fA[0] = -2539; // [cm]
227 fA[1] = 900; // [cm]
228 fA[2] = 17.5; // [cm]
229 fA[3] = 1;// 4;
230 fA[4] = 58.3;
231
232 TArrayD fB(7);
233 fB[0] = -8.64; // [GeV]
234 fB[1] = -0.069; // [GeV]
235 fB[2] = 0.00066; // [GeV]
236 fB[3] = 0.033; // [GeV]
237 fB[4] = 0.000226; // [GeV]
238 fB[5] = 4.14e-8; // [GeV]
239 fB[6] = -0.06;
240
241 // Set starting values and step sizes for parameters
242 for (Int_t i=0; i<fA.GetSize(); i++)
243 {
244 TString name = "fA[";
245 name += i;
246 name += "]";
247 Double_t vinit = fA[i];
248 Double_t step = fabs(fA[i]/3);
249
250 Double_t limlo = 0; // limlo=limup=0: no limits
251 Double_t limup = 0;
252
253 Bool_t rc = minuit.DefineParameter(i, name, vinit, step, limlo, limup);
254 if (evalenergy)
255 minuit.FixParameter(i);
256
257 if (!rc)
258 continue;
259
260 cout << "Error in defining parameter #" << i << endl;
261 return;
262 }
263
264 for (Int_t i=0; i<fB.GetSize(); i++)
265 {
266 TString name = "fB[";
267 name += i;
268 name += "]";
269 Double_t vinit = fB[i];
270 Double_t step = fabs(fB[i]/3);
271
272 Double_t limlo = 0; // limlo=limup=0: no limits
273 Double_t limup = 0;
274
275 Bool_t rc = minuit.DefineParameter(i+fA.GetSize(), name, vinit, step, limlo, limup);
276 if (!evalenergy)
277 minuit.FixParameter(i+fA.GetSize());
278
279 if (!rc)
280 continue;
281
282 cout << "Error in defining parameter #" << i+fA.GetSize() << endl;
283 return;
284 }
285
286 //
287 // Setup globals used in FCN
288 //
289 minuit.SetObjectFit(&evtloop);
290
291 TStopwatch clock;
292 clock.Start();
293
294 // Now ready for minimization step: minuit.mnexcm("MIGRAD", arglist, 1, ierflg)
295 gLog.SetNullOutput(kTRUE);
296 Bool_t rc = minuit.Migrad();
297 gLog.SetNullOutput(kFALSE);
298
299 if (rc)
300 {
301 cout << "Migrad failed." << endl;
302 return;
303 }
304
305 cout << endl;
306 clock.Stop();
307 clock.Print();
308 cout << endl;
309
310 for (Int_t i=(evalenergy?fA.GetSize():0); i<(evalenergy?fA.GetSize()+fB.GetSize():fA.GetSize()); i++)
311 {
312 Double_t val;
313 Double_t er;
314
315 if (!minuit.GetParameter(i, val, er))
316 {
317 cout << "Error getting parameter #" << i << endl;
318 return;
319 }
320
321 cout << i << ": " << val << " +- " << er << endl;
322 }
323
324 /*
325 // Print results
326 Double_t amin, edm, errdef;
327 Int_t nvpar, nparx, icstat;
328 minuit.mnstat(amin, edm, errdef, nvpar, nparx, icstat);
329 minuit.mnprin(3, amin);
330 */
331}
Note: See TracBrowser for help on using the repository browser.