source: trunk/MagicSoft/Mars/macros/CT1EnergyEst.C@ 1852

Last change on this file since 1852 was 1852, checked in by moralejo, 22 years ago
*** empty log message ***
File size: 10.4 KB
Line 
1/* ======================================================================== *\
2!
3! *
4! * This file is part of MARS, the MAGIC Analysis and Reconstruction
5! * Software. It is distributed to you in the hope that it can be a useful
6! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
7! * It is distributed WITHOUT ANY WARRANTY.
8! *
9! * Permission to use, copy, modify and distribute this software and its
10! * documentation for any purpose is hereby granted without fee,
11! * provided that the above copyright notice appear in all copies and
12! * that both that copyright notice and this permission notice appear
13! * in supporting documentation. It is provided "as is" without express
14! * or implied warranty.
15! *
16!
17!
18! Author(s): Thomas Bretz et al, 09/2002 <mailto:tbretz@astro.uni-wuerzburg.de>
19!
20! Copyright: MAGIC Software Development, 2000-2002
21!
22!
23\* ======================================================================== */
24
25void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
26{
27 MEvtLoop *evtloop = (MEvtLoop*)gMinuit->GetObjectFit();
28
29 MTaskList *tlist = (MTaskList*)evtloop->GetParList()->FindObject("MTaskList"); //GetTaskList();
30
31 MChisqEval *eval = (MChisqEval*) tlist->FindObject("MChisqEval");
32 MEnergyEstParam *eest = (MEnergyEstParam*)tlist->FindObject("MEnergyEstParam");
33
34 eest->SetCoeff(TArrayD(eest->GetNumCoeff(), par));
35
36 evtloop->Eventloop();
37
38 f = eval->GetChisq();
39}
40
41//------------------------------------------------------------------------
42
43void CT1EnergyEst(Float_t maxhadronness=1.)
44{
45 // Bool_t evalenergy=kFALSE;
46 //
47 // Fill events into a MHMatrix
48 //
49 MParList parlist;
50 MHMatrix matrix;
51
52 Int_t col = matrix.AddColumn("MMcEvt.fEnergy");
53
54 MEnergyEstParam eest("Hillas");
55 eest.Add("HillasSrc");
56 eest.InitMapping(&matrix);
57
58 MReadTree read("Events");
59 read.AddFile("MC_ON2.root", 200);
60 read.DisableAutoScheme();
61
62 TString hcut("MHadronness.fHadronness<");
63 hcut += maxhadronness;
64
65 MF filterhadrons(hcut);
66
67 if (!matrix.Fill(&parlist, &read))
68 return;
69
70 //
71 // Setup the tasklist used to evaluate the needed chisq
72 //
73 MTaskList tasklist;
74 parlist.AddToList(&tasklist);
75
76 MMatrixLoop loop(&matrix);
77
78
79 MChisqEval eval;
80 eval.SetY1(new MDataElement(&matrix, col));
81 eval.SetY2(new MDataMember("MEnergyEst.fEnergy"));
82 eval.SetOwner();
83
84 tasklist.AddToList(&loop);
85 tasklist.AddToList(&eest);
86 tasklist.AddToList(&eval);
87
88 MEvtLoop evtloop;
89 evtloop.SetParList(&parlist);
90
91 //
92 // Be careful: This is not thread safe
93 //
94 TMinuit minuit(12);
95 minuit.SetPrintLevel(-1);
96 minuit.SetFCN(fcn);
97
98 // Ready for: minuit.mnexcm("SET ERR", arglist, 1, ierflg)
99 if (minuit.SetErrorDef(1))
100 {
101 cout << "SetErrorDef failed." << endl;
102 return;
103 }
104
105 //
106 // Set initial values
107 //
108 TArrayD fA(5);
109 fA[0] =1;//4916.4; */-2414.75;
110 fA[1] =1;//149.549; */ 1134.28;
111 fA[2] =1;//-558.209; */ 132.932;
112 fA[3] =1;//0.270725; */ 0.292845;
113 fA[4] =1;//107.001; */ 107.001;
114
115 TArrayD fB(7);
116 fB[0] = 1;//-8234.12; */ -4282.25;
117 fB[1] = 1;// 23.2153; */ 18.892;
118 fB[2] = 1;// 0.416372;*/ 0.193373;
119 fB[3] = 1;// 332.42; */ 203.803;
120 fB[4] = 1;// -0.701764;*/ -0.534876;
121 fB[5] = 1;//-0.0131774;*/ -0.00789539;
122 fB[6] = 1;//-0.162687;*/ 0.111913;
123
124 // Set starting values and step sizes for parameters
125 for (Int_t i=0; i<fA.GetSize(); i++)
126 {
127 TString name = "fA[";
128 name += i;
129 name += "]";
130 Double_t vinit = fA[i];
131 Double_t step = fabs(fA[i]/3);
132
133 Double_t limlo = 0; // limlo=limup=0: no limits
134 Double_t limup = 0;
135
136 Bool_t rc = minuit.DefineParameter(i, name, vinit, step, limlo, limup);
137 if (!rc)
138 continue;
139
140 cout << "Error in defining parameter #" << i << endl;
141 return;
142 }
143
144 for (Int_t i=0; i<fB.GetSize(); i++)
145 {
146 TString name = "fB[";
147 name += i;
148 name += "]";
149 Double_t vinit = fB[i];
150 Double_t step = fabs(fB[i]/3);
151
152 Double_t limlo = 0; // limlo=limup=0: no limits
153 Double_t limup = 0;
154
155 Bool_t rc = minuit.DefineParameter(i+fA.GetSize(), name, vinit, step, limlo, limup);
156 if (!rc)
157 continue;
158
159 cout << "Error in defining parameter #" << i+fA.GetSize() << endl;
160 return;
161 }
162
163 //
164 // Setup globals used in FCN
165 //
166 minuit.SetObjectFit(&evtloop);
167
168
169 TStopwatch clock;
170 clock.Start();
171
172 // Now ready for minimization step: minuit.mnexcm("MIGRAD", arglist, 1, ierflg)
173 gLog.SetNullOutput(kTRUE);
174 Bool_t rc = minuit.Migrad();
175 gLog.SetNullOutput(kFALSE);
176
177 if (rc)
178 {
179 cout << "Migrad failed." << endl;
180 return;
181 }
182
183 cout << endl;
184 clock.Stop();
185 clock.Print();
186 cout << endl;
187
188 cout << "Resulting Chisq: " << minuit.fAmin << endl;
189
190 for (Int_t i=0; i<fA.GetSize()+fB.GetSize(); i++)
191 {
192 Double_t val;
193 Double_t er;
194
195 if (!minuit.GetParameter(i, val, er))
196 {
197 cout << "Error getting parameter #" << i << endl;
198 return;
199 }
200
201 cout << i << ": " << val << " +- " << er << endl;
202 }
203
204 /*
205 // Print results
206 Double_t amin, edm, errdef;
207 Int_t nvpar, nparx, icstat;
208 minuit.mnstat(amin, edm, errdef, nvpar, nparx, icstat);
209 minuit.mnprin(3, amin);
210 */
211 eest.Print();
212
213
214
215 //---------------------------------------------------------------
216
217 // Part 2: Now test how the energy estimation method works.
218 //
219 //
220 // Create a empty Parameter List and an empty Task List
221 // The tasklist is identified in the eventloop by its name
222 //
223
224
225 MTaskList tlist2;
226 parlist.Replace(&tlist2);
227
228 //
229 // Now setup the tasks and tasklist:
230 // ---------------------------------
231 //
232
233
234 read->SetEventNum(0);
235
236 //
237 // Use this to change the binnign of the histograms to CT1-style
238 //
239 Bool_t usect1 = kTRUE;
240
241 MH3 mh3e("MMcEvt.fEnergy", "(MEnergyEst.fEnergy/MMcEvt.fEnergy-1)*(MEnergyEst.fEnergy/MMcEvt.fEnergy-1)");
242 MH3 mh3i("MMcEvt.fImpact/100", "(MEnergyEst.fImpact/MMcEvt.fImpact-1)*(MEnergyEst.fImpact/MMcEvt.fImpact-1)");
243 MH3 mh3eo("MMcEvt.fEnergy", "MEnergyEst.fEnergy/MMcEvt.fEnergy-1");
244 MH3 mh3io("MMcEvt.fImpact/100", "MEnergyEst.fImpact/MMcEvt.fImpact-1");
245
246 MH3 mh3e2("MEnergyEst.fEnergy", "(MEnergyEst.fEnergy/MMcEvt.fEnergy-1)*(MEnergyEst.fEnergy/MMcEvt.fEnergy-1)");
247 MH3 mh3i2("MEnergyEst.fImpact/100", "(MEnergyEst.fImpact/MMcEvt.fImpact-1)*(MEnergyEst.fImpact/MMcEvt.fImpact-1)");
248 MH3 mh3eo2("MEnergyEst.fEnergy", "MEnergyEst.fEnergy/MMcEvt.fEnergy-1");
249 MH3 mh3io2("MEnergyEst.fImpact/100", "MEnergyEst.fImpact/MMcEvt.fImpact-1");
250
251 MH3 mhe("MMcEvt.fEnergy", "MEnergyEst.fEnergy");
252 MH3 mhi("MMcEvt.fImpact/100", "MEnergyEst.fImpact/100");
253
254 mh3e.SetName("HistEnergy");
255 mh3i.SetName("HistImpact");
256 mh3eo.SetName("HistEnergyOffset");
257 mh3io.SetName("HistImpactOffset");
258
259 mh3e2.SetName("HistEnergy");
260 mh3i2.SetName("HistImpact");
261 mh3eo2.SetName("HistEnergyOffset");
262 mh3io2.SetName("HistImpactOffset");
263
264 mhe.SetName("HistEE");
265 mhi.SetName("HistII");
266
267 MFillH hfille(&mh3e);
268 MFillH hfilli(&mh3i);
269 MFillH hfilleo(&mh3eo);
270 MFillH hfillio(&mh3io);
271
272 MFillH hfille2(&mh3e2);
273 MFillH hfilli2(&mh3i2);
274 MFillH hfilleo2(&mh3eo2);
275 MFillH hfillio2(&mh3io2);
276
277 MFillH hfillee(&mhe);
278 MFillH hfillii(&mhi);
279
280 MBinning binsex("BinningHistEnergyX");
281 MBinning binsey("BinningHistEnergyY");
282 MBinning binsix("BinningHistImpactX");
283 MBinning binsiy("BinningHistImpactY");
284 MBinning binseox("BinningHistEnergyOffsetX");
285 MBinning binseoy("BinningHistEnergyOffsetY");
286 MBinning binsiox("BinningHistImpactOffsetX");
287 MBinning binsioy("BinningHistImpactOffsetY");
288 MBinning binseex("BinningHistEEX");
289 MBinning binsiix("BinningHistIIX");
290 MBinning binseey("BinningHistEEY");
291 MBinning binsiiy("BinningHistIIY");
292
293 binsex.SetEdgesLog(50, usect1 ? 300: 10, usect1 ? 50000 : 1e4);
294 binsey.SetEdges(50, 0, usect1 ? 0.8 : 1.75);
295 binseox.SetEdgesLog(50, usect1 ? 300 : 10, usect1 ? 50000 : 1e4);
296 binseoy.SetEdges(50, usect1 ? -0.75 : -1.75, usect1 ? 0.75 : 1.75);
297
298 binsix.SetEdges(50, 0, usect1 ? 275 : 300);
299 binsiy.SetEdges(50, 0, usect1 ? 0.2 : 1.75);
300 binsiox.SetEdges(50, 0, usect1 ? 275 : 300);
301 binsioy.SetEdges(50, usect1 ? -0.75 : -1.75, usect1 ? 0.75 : 1.75);
302
303 binseex.SetEdgesLog(50, usect1 ? 300 : 10, usect1 ? 50000 : 15e3);
304 binseey.SetEdgesLog(50, usect1 ? 300 : 1, usect1 ? 50000 : 2e3);
305 binsiix.SetEdges(50, 0, usect1 ? 275 : 300);
306 binsiiy.SetEdges(50, 0, usect1 ? 275 : 150);
307
308 parlist.AddToList(&binsex);
309 parlist.AddToList(&binsey);
310 parlist.AddToList(&binsix);
311 parlist.AddToList(&binsiy);
312 parlist.AddToList(&binseox);
313 parlist.AddToList(&binseoy);
314 parlist.AddToList(&binsiox);
315 parlist.AddToList(&binsioy);
316 parlist.AddToList(&binseex);
317 parlist.AddToList(&binseey);
318 parlist.AddToList(&binsiix);
319 parlist.AddToList(&binsiiy);
320
321 eest.StopMapping();
322 eest.Add("HillasSrc");
323
324 //
325 // Setup tasklists
326 //
327 tlist2.AddToList(&read);
328
329 TString hcut2("MHadronness.fHadronness>");
330 hcut2 += maxhadronness;
331 MContinue cont(hcut2);
332
333 tlist2.AddToList(&cont);
334
335 tlist2.AddToList(&eest);
336 // tlist2.AddToList(new MPrint("MMcEvt"));
337 // tlist2.AddToList(new MPrint("MEnergyEst"));
338
339 tlist2.AddToList(&hfille);
340 tlist2.AddToList(&hfilli);
341 tlist2.AddToList(&hfilleo);
342 tlist2.AddToList(&hfillio);
343
344 tlist2.AddToList(&hfille2);
345 tlist2.AddToList(&hfilli2);
346 tlist2.AddToList(&hfilleo2);
347 tlist2.AddToList(&hfillio2);
348
349 tlist2.AddToList(&hfillee);
350 tlist2.AddToList(&hfillii);
351
352 //
353 // Create and setup the eventloop
354 //
355 MProgressBar bar;
356
357 MEvtLoop evtloop2;
358 evtloop2.SetProgressBar(&bar);
359 evtloop2.SetParList(&parlist);
360
361 //
362 // Execute your analysis
363 //
364 if (!evtloop2.Eventloop())
365 return;
366
367 tlist2.PrintStatistics();
368
369 const TString text = "\\sqrt{<y>}=%.0f%%";
370
371 char txt[1000];
372
373 TCanvas *c=new TCanvas("Est1", "Estimates vs. E_{true}");
374 c->Divide(2,2);
375 c->cd(1);
376 mh3i.DrawClone("PROFXnonew");
377 sprintf(txt, text.Data(), sqrt(mh3i.GetHist().GetMean(2))*100);
378 TLatex *t = new TLatex(180, 0.15, txt);
379 t->Draw();
380 c->cd(2);
381 mh3e.DrawClone("PROFXnonew");
382 sprintf(txt, text.Data(), sqrt(mh3e.GetHist().GetMean(2))*100);
383 t = new TLatex(3.5, 0.6, txt);
384 t->Draw();
385 c->cd(3);
386 mh3io.DrawClone("PROFXnonew");
387 c->cd(4);
388 mh3eo.DrawClone("PROFXnonew");
389
390 c=new TCanvas("Est2", "Estimates vs. E_{est}");
391 c->Divide(2,2);
392 c->cd(1);
393 mh3i2.DrawClone("PROFXnonew");
394 c->cd(2);
395 mh3e2.DrawClone("PROFXnonew");
396 c->cd(3);
397 mh3io2.DrawClone("PROFXnonew");
398 c->cd(4);
399 mh3eo2.DrawClone("PROFXnonew");
400
401 mhe.DrawClone("PROFX");
402 mhi.DrawClone("PROFX");
403}
Note: See TracBrowser for help on using the repository browser.