source: trunk/MagicSoft/Mars/mjtrain/MJTrainDisp.cc@ 8672

Last change on this file since 8672 was 8657, checked in by tbretz, 17 years ago
*** empty log message ***
File size: 14.2 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 11/2005 <mailto:tbretz@astro.uni-wuerzburg.de>
19!
20! Copyright: MAGIC Software Development, 2005-2007
21!
22!
23\* ======================================================================== */
24
25/////////////////////////////////////////////////////////////////////////////
26//
27// MJTrainDisp
28//
29//
30// Example:
31// --------
32//
33// // SequencesOn are used for training, SequencesOff for Testing
34// MDataSet set("mctesttrain.txt");
35// set.SetNumAnalysis(1); // Must have a number
36// MJTrainDisp opt;
37// //opt.SetDebug();
38// opt.AddParameter("MHillas.fLength");
39// opt.AddParameter("MHillas.fWidth");
40// MStatusDisplay *d = new MStatusDisplay;
41// opt.SetDisplay(d);
42// opt.AddPreCut("MHillasSrc.fDCA*MGeomCam.fConvMm2Deg<0.3");
43// opt.Train("rf-disp.root", set, 30000); // Number of train events
44//
45//
46// Random Numbers:
47// ---------------
48// Use:
49// if(gRandom)
50// delete gRandom;
51// gRandom = new TRandom3();
52// in advance to change the random number generator.
53//
54////////////////////////////////////////////////////////////////////////////
55#include "MJTrainDisp.h"
56
57#include <TLine.h>
58#include <TCanvas.h>
59
60#include "MHMatrix.h"
61
62#include "MLog.h"
63#include "MLogManip.h"
64
65// tools
66#include "MDataSet.h"
67#include "MTFillMatrix.h"
68#include "MChisqEval.h"
69#include "MStatusDisplay.h"
70
71// eventloop
72#include "MParList.h"
73#include "MTaskList.h"
74#include "MEvtLoop.h"
75
76// tasks
77#include "MReadMarsFile.h"
78#include "MContinue.h"
79#include "MFillH.h"
80#include "MRanForestCalc.h"
81#include "MParameterCalc.h"
82
83// container
84#include "MParameters.h"
85
86// histograms
87#include "MBinning.h"
88#include "MH3.h"
89#include "MHThetaSq.h"
90
91// filter
92#include "MFilterList.h"
93
94ClassImp(MJTrainDisp);
95
96using namespace std;
97
98const TString MJTrainDisp::fgTrainParameter = "MHillasSrc.fDist*MGeomCam.fConvMm2Deg";
99
100// --------------------------------------------------------------------------
101//
102// Display a result histogram either vs. size or energy
103// FIXME: Could be moved into a new histogram class.
104//
105void MJTrainDisp::DisplayHist(TCanvas &c, Int_t i, MH3 &mh3) const
106{
107 MH::SetPalette("pretty");
108
109 TH1 &hist = *(TH1*)mh3.GetHist().Clone();
110 hist.SetBit(TH1::kNoStats);
111 hist.SetDirectory(0);
112
113 TLine line;
114 line.SetLineStyle(kDashed);
115 line.SetLineWidth(1);
116
117 c.cd(i);
118 gPad->SetBorderMode(0);
119 gPad->SetFrameBorderMode(0);
120 //gPad->SetFillColor(kWhite);
121 gPad->SetLogx();
122 gPad->SetGridx();
123 gPad->SetGridy();
124 //gPad->SetLeftMargin(0.12);
125 //gPad->SetRightMargin(0.12);
126
127 for (int x=0; x<=hist.GetNbinsX(); x++)
128 {
129 Float_t n = 0;
130 for (int y=1; y<=2; y++)
131 n += hist.GetBinContent(x,y);
132
133 if (n==0)
134 continue;
135
136 for (int y=0; y<=hist.GetNbinsY(); y++)
137 hist.SetBinContent(x, y, 200*hist.GetBinContent(x,y)/n);
138 }
139
140 hist.SetMaximum(100);
141 hist.DrawCopy("colz");
142
143 line.DrawLine(10, 0.04, 31623, 0.04);
144
145 c.cd(i+2);
146 gPad->SetBorderMode(0);
147 gPad->SetFrameBorderMode(0);
148 //gPad->SetFillColor(kWhite);
149 gPad->SetLogx();
150 gPad->SetGridx();
151 gPad->SetGridy();
152 //gPad->SetLeftMargin(0.12);
153 //gPad->SetRightMargin(0.12);
154
155 for (int x=0; x<=hist.GetNbinsX(); x++)
156 {
157 Float_t n = 0;
158 for (int y=0; y<=hist.GetNbinsY(); y++)
159 n += hist.GetBinContent(x,y);
160
161 if (n==0)
162 continue;
163
164 for (int y=0; y<=hist.GetNbinsY(); y++)
165 hist.SetBinContent(x, y, 100*hist.GetBinContent(x,y)/n);
166 }
167
168 hist.SetMaximum(25);
169 hist.DrawCopy("colz");
170
171 line.DrawLine(10, 0.04, 31623, 0.04);
172}
173
174// --------------------------------------------------------------------------
175//
176// Display the result histograms in a new tab.
177//
178void MJTrainDisp::DisplayResult(MH3 &hsize, MH3 &henergy)
179{
180 TCanvas &c = fDisplay->AddTab("Disp");
181 c.Divide(2,2);
182
183 DisplayHist(c, 1, hsize);
184 DisplayHist(c, 2, henergy);
185}
186
187// --------------------------------------------------------------------------
188//
189// Run Disp optimization
190//
191Bool_t MJTrainDisp::Train(const char *out, const MDataSet &set, Int_t num)
192{
193 SetTitle(Form("TrainDisp: %s", out));
194
195 if (fDisplay)
196 fDisplay->SetTitle(fTitle);
197
198 if (!set.IsValid())
199 {
200 *fLog << err << "ERROR - DataSet invalid!" << endl;
201 return kFALSE;
202 }
203
204 if (!HasWritePermission(out))
205 return kFALSE;
206
207 *fLog << inf;
208 fLog->Separator(GetDescriptor());
209
210 // --------------------- Setup files --------------------
211 MReadMarsFile readtrn("Events");
212 MReadMarsFile readtst("Events");
213 readtrn.DisableAutoScheme();
214 readtst.DisableAutoScheme();
215
216 if (!set.AddFilesOn(readtrn))
217 return kFALSE;
218 if (!set.AddFilesOff(readtst))
219 return kFALSE;
220
221 // ----------------------- Setup Matrix ------------------
222 MHMatrix train("Train");
223 train.AddColumns(fRules);
224 if (fEnableWeights)
225 train.AddColumn("MWeight.fVal");
226 train.AddColumn(fTrainParameter);
227
228 // ----------------------- Fill Matrix RF ----------------------
229 MTFillMatrix fill(fTitle);
230 fill.SetDisplay(fDisplay);
231 fill.SetLogStream(fLog);
232 fill.SetDestMatrix1(&train, num);
233 fill.SetReader(&readtrn);
234 fill.AddPreCuts(fPreCuts);
235 fill.AddPreCuts(fTrainCuts);
236 fill.AddPreTasks(fPreTasks);
237 fill.AddPostTasks(fPostTasks);
238 if (!fill.Process())
239 return kFALSE;
240
241 // ------------------------ Train RF --------------------------
242 MRanForestCalc rf("TrainDisp", fTitle);
243 rf.SetNumTrees(fNumTrees);
244 rf.SetNdSize(fNdSize);
245 rf.SetNumTry(fNumTry);
246 rf.SetNumObsoleteVariables(1);
247 rf.SetLastDataColumnHasWeights(fEnableWeights);
248 rf.SetDisplay(fDisplay);
249 rf.SetLogStream(fLog);
250 rf.SetFileName(out);
251 rf.SetDebug(fDebug>1);
252 rf.SetNameOutput("Disp");
253
254 /*
255 MBinning b(32, 10, 100000, "BinningEnergyEst", "log");
256
257 if (!rf.TrainMultiRF(train, b.GetEdgesD())) // classification with one tree per bin
258 return;
259
260 if (!rf.TrainSingleRF(train, b.GetEdgesD())) // classification into different bins
261 return;
262 */
263 if (!rf.TrainRegression(train)) // regression (best choice)
264 return kFALSE;
265
266 // --------------------- Display result ----------------------
267
268 gLog.Separator("Test");
269
270 MH::SetPalette("pretty");
271
272 MParList plist;
273 MTaskList tlist;
274 plist.AddToList(this);
275 plist.AddToList(&tlist);
276 //plist.AddToList(&b);
277
278 MParameterD par("ThetaSquaredCut");
279 par.SetVal(0.2);
280 plist.AddToList(&par);
281
282 MAlphaFitter fit;
283 fit.SetPolynomOrder(0);
284 fit.SetSignalFitMax(0.8);
285 fit.EnableBackgroundFit(kFALSE);
286 fit.SetSignalFunction(MAlphaFitter::kThetaSq);
287 fit.SetMinimizationStrategy(MAlphaFitter::kGaussSigma);
288 plist.AddToList(&fit);
289
290 MFilterList list;
291 if (!list.AddToList(fPreCuts))
292 *fLog << err << "ERROR - Calling MFilterList::AddToList for fPreCuts failed!" << endl;
293 if (!list.AddToList(fTestCuts))
294 *fLog << err << "ERROR - Calling MFilterList::AddToList for fTestCuts failed!" << endl;
295
296 MContinue cont(&list);
297 cont.SetInverted();
298
299 MHThetaSq hist;
300 hist.SkipHistTime();
301 hist.SkipHistTheta();
302 hist.SkipHistEnergy();
303
304 MFillH fillh(&hist, "", "FillThetaSq");
305
306 // 0 = disp^2 - 2*disp*dist*cos(alpha) + dist^2
307
308 // cos^2 -1 = - sin^2
309
310 // disp = +dist* (cos(alpha) +/- sqrt(cos^2(alpha) - 1) )
311
312 const char *rule = "(MHillasSrc.fDist*MGeomCam.fConvMm2Deg)^2 + (Disp.fVal)^2 - (2*MHillasSrc.fDist*MGeomCam.fConvMm2Deg*Disp.fVal*cos(MHillasSrc.fAlpha*kDeg2Rad))";
313
314 MParameterCalc calcthetasq(rule, "MThetaSqCalc");
315 calcthetasq.SetNameParameter("ThetaSquared");
316
317 MChisqEval eval;
318 eval.SetY1("sqrt(ThetaSquared.fVal)");
319
320 MH3 hdisp1("MHillas.fSize", "ThetaSquared.fVal");
321 MH3 hdisp2("MMcEvt.fEnergy", "ThetaSquared.fVal");
322 hdisp1.SetTitle("\\vartheta distribution vs. Size:Size [phe]:\\vartheta [\\circ]");
323 hdisp2.SetTitle("\\vartheta distribution vs. Energy:Enerhy [GeV]:\\vartheta [\\circ]");
324
325 MBinning binsx( 70, 10, 31623, "BinningMH3X", "log");
326 MBinning binsy( 75, 0, 0.3, "BinningMH3Y", "lin");
327
328 plist.AddToList(&binsx);
329 plist.AddToList(&binsy);
330
331 MFillH fillh2a(&hdisp1, "", "FillSize");
332 MFillH fillh2b(&hdisp2, "", "FillEnergy");
333 fillh2a.SetDrawOption("colz profx");
334 fillh2b.SetDrawOption("colz profx");
335 fillh2a.SetNameTab("Size");
336 fillh2b.SetNameTab("Energy");
337
338 tlist.AddToList(&readtst);
339 tlist.AddToList(&cont);
340 tlist.AddToList(&rf);
341 tlist.AddToList(&calcthetasq);
342 tlist.AddToList(&fillh);
343 tlist.AddToList(&fillh2a);
344 tlist.AddToList(&fillh2b);
345 tlist.AddToList(&eval);
346
347 MEvtLoop loop(fTitle);
348 loop.SetLogStream(fLog);
349 loop.SetDisplay(fDisplay);
350 loop.SetParList(&plist);
351 //if (!SetupEnv(loop))
352 // return kFALSE;
353
354 if (!loop.Eventloop())
355 return kFALSE;
356
357 // Print the result
358 *fLog << inf;
359 *fLog << "Rule: " << rule << endl;
360 *fLog << "Disp: " << fTrainParameter << endl;
361 hist.GetAlphaFitter().Print("result");
362
363 DisplayResult(hdisp1, hdisp2);
364
365 SetPathOut(out);
366 if (!WriteDisplay(0, "UPDATE"))
367 return kFALSE;
368
369 return kTRUE;
370}
371
372/*
373#include "MParameterCalc.h"
374#include "MHillasCalc.h"
375#include "../mpointing/MSrcPosRndm.h"
376
377Bool_t MJTrainDisp::TrainGhostbuster(const char *out, const MDataSet &set, Int_t num)
378{
379 SetTitle(Form("TrainGhostbuster: %s", out));
380
381 if (fDisplay)
382 fDisplay->SetTitle(fTitle);
383
384 if (!set.IsValid())
385 {
386 *fLog << err << "ERROR - DataSet invalid!" << endl;
387 return kFALSE;
388 }
389
390 if (!HasWritePermission(out))
391 return kFALSE;
392
393 *fLog << inf;
394 fLog->Separator(GetDescriptor());
395
396 // --------------------- Setup files --------------------
397 MReadMarsFile readtrn("Events");
398 MReadMarsFile readtst("Events");
399 readtrn.DisableAutoScheme();
400 readtst.DisableAutoScheme();
401
402 if (!set.AddFilesOn(readtrn))
403 return kFALSE;
404 if (!set.AddFilesOff(readtst))
405 return kFALSE;
406
407 // ----------------------- Setup Matrix ------------------
408 MHMatrix train("Train");
409 train.AddColumns(fRules);
410 if (fEnableWeights)
411 train.AddColumn("MWeight.fVal");
412 train.AddColumn("sign(MHillasSrc.fCosDeltaAlpha)==sign(SignStore.fVal)");
413
414 MParameterCalc calc("MHillasSrc.fCosDeltaAlpha", "SignStore");
415 calc.SetNameParameter("SignStore");
416
417 MSrcPosRndm rndm;
418 rndm.SetRule(fTrainParameter);
419 //rndm.SetDistOfSource(120*3.37e-3);
420
421 MHillasCalc hcalc;
422 hcalc.SetFlags(MHillasCalc::kCalcHillasSrc);
423
424 // ----------------------- Fill Matrix RF ----------------------
425 MTFillMatrix fill(fTitle);
426 fill.SetDisplay(fDisplay);
427 fill.SetLogStream(fLog);
428 fill.SetDestMatrix1(&train, num);
429 fill.SetReader(&readtrn);
430 fill.AddPreCuts(fPreCuts);
431 fill.AddPreCuts(fTrainCuts);
432 fill.AddPreTasks(fPreTasks);
433 fill.AddPostTasks(fPostTasks);
434 fill.AddPostTask(&calc);
435 fill.AddPostTask(&rndm);
436 fill.AddPostTask(&hcalc);
437 if (!fill.Process())
438 return kFALSE;
439
440 // ------------------------ Train RF --------------------------
441 MRanForestCalc rf("TrainGhostbuster", fTitle);
442 rf.SetNumTrees(fNumTrees);
443 rf.SetNdSize(fNdSize);
444 rf.SetNumTry(fNumTry);
445 rf.SetNumObsoleteVariables(1);
446 rf.SetLastDataColumnHasWeights(fEnableWeights);
447 rf.SetDisplay(fDisplay);
448 rf.SetLogStream(fLog);
449 rf.SetFileName(out);
450 rf.SetDebug(fDebug>1);
451 rf.SetNameOutput("Sign");
452
453 if (!rf.TrainRegression(train)) // regression (best choice)
454 return kFALSE;
455
456 // --------------------- Display result ----------------------
457
458 gLog.Separator("Test");
459
460 MH::SetPalette("pretty");
461
462 MParList plist;
463 MTaskList tlist;
464 plist.AddToList(this);
465 plist.AddToList(&tlist);
466 //plist.AddToList(&b);
467
468 MFilterList list;
469 if (!list.AddToList(fPreCuts))
470 *fLog << err << "ERROR - Calling MFilterList::AddToList for fPreCuts failed!" << endl;
471 if (!list.AddToList(fTestCuts))
472 *fLog << err << "ERROR - Calling MFilterList::AddToList for fTestCuts failed!" << endl;
473
474 MContinue cont(&list);
475 cont.SetInverted();
476
477 const char *rule = "abs(Sign.fVal-(sign(MHillasSrc.fCosDeltaAlpha)==sign(SignStore.fVal)))";
478
479 //MChisqEval eval;
480 //eval.SetY1("sqrt(ThetaSquared.fVal)");
481
482 MH3 hsign1("MHillas.fSize", rule);
483 MH3 hsign2("MMcEvt.fEnergy", rule);
484 hsign1.SetTitle("Was ist das? vs. Size:Size [phe]:XXX");
485 hsign2.SetTitle("Was ist das? vs. Energy:Enerhy [GeV]:XXXX");
486
487 MBinning binsx( 70, 10, 31623, "BinningMH3X", "log");
488 MBinning binsy( 51, 0, 1, "BinningMH3Y", "lin");
489
490 plist.AddToList(&binsx);
491 plist.AddToList(&binsy);
492
493 MFillH fillh2a(&hsign1, "", "FillSize");
494 MFillH fillh2b(&hsign2, "", "FillEnergy");
495 fillh2a.SetDrawOption("profx colz");
496 fillh2b.SetDrawOption("profx colz");
497 fillh2a.SetNameTab("Size");
498 fillh2b.SetNameTab("Energy");
499
500 tlist.AddToList(&readtst);
501 tlist.AddToList(&cont);
502 tlist.AddToList(&calc);
503 tlist.AddToList(&rndm);
504 tlist.AddToList(&hcalc);
505 tlist.AddToList(&rf);
506 tlist.AddToList(&fillh2a);
507 tlist.AddToList(&fillh2b);
508 //tlist.AddToList(&eval);
509
510 MEvtLoop loop(fTitle);
511 loop.SetLogStream(fLog);
512 loop.SetDisplay(fDisplay);
513 loop.SetParList(&plist);
514 //if (!SetupEnv(loop))
515 // return kFALSE;
516
517 if (!loop.Eventloop())
518 return kFALSE;
519
520 // Print the result
521 *fLog << inf << "Rule: " << rule << endl;
522
523 //DisplayResult(hdisp1, hdisp2);
524
525 SetPathOut(out);
526 if (!WriteDisplay(0, "UPDATE"))
527 return kFALSE;
528
529 return kTRUE;
530}
531*/
Note: See TracBrowser for help on using the repository browser.