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 | // // Two of the displayed histograms show the cut efficiency for
|
---|
46 | // // a given Theta-Cut. The default is 0.215. It can be overwritten
|
---|
47 | // // by
|
---|
48 | // opt.SetThetaCut(0.165);
|
---|
49 | //
|
---|
50 | //
|
---|
51 | // Random Numbers:
|
---|
52 | // ---------------
|
---|
53 | // Use:
|
---|
54 | // if(gRandom)
|
---|
55 | // delete gRandom;
|
---|
56 | // gRandom = new TRandom3();
|
---|
57 | // in advance to change the random number generator.
|
---|
58 | //
|
---|
59 | ////////////////////////////////////////////////////////////////////////////
|
---|
60 | #include "MJTrainDisp.h"
|
---|
61 |
|
---|
62 | #include <TLine.h>
|
---|
63 | #include <TCanvas.h>
|
---|
64 |
|
---|
65 | #include "MHMatrix.h"
|
---|
66 |
|
---|
67 | #include "MLog.h"
|
---|
68 | #include "MLogManip.h"
|
---|
69 |
|
---|
70 | // tools
|
---|
71 | #include "MDataSet.h"
|
---|
72 | #include "MTFillMatrix.h"
|
---|
73 | #include "MChisqEval.h"
|
---|
74 | #include "MStatusDisplay.h"
|
---|
75 |
|
---|
76 | // eventloop
|
---|
77 | #include "MParList.h"
|
---|
78 | #include "MTaskList.h"
|
---|
79 | #include "MEvtLoop.h"
|
---|
80 |
|
---|
81 | // tasks
|
---|
82 | #include "MReadMarsFile.h"
|
---|
83 | #include "MContinue.h"
|
---|
84 | #include "MFillH.h"
|
---|
85 | #include "MRanForestCalc.h"
|
---|
86 | #include "MParameterCalc.h"
|
---|
87 |
|
---|
88 | // container
|
---|
89 | #include "MParameters.h"
|
---|
90 |
|
---|
91 | // histograms
|
---|
92 | #include "MBinning.h"
|
---|
93 | #include "MH3.h"
|
---|
94 | #include "MHn.h"
|
---|
95 | #include "MHThetaSq.h"
|
---|
96 |
|
---|
97 | // filter
|
---|
98 | #include "MFilterList.h"
|
---|
99 |
|
---|
100 | ClassImp(MJTrainDisp);
|
---|
101 |
|
---|
102 | using namespace std;
|
---|
103 |
|
---|
104 | const TString MJTrainDisp::fgTrainParameter = "MHillasSrc.fDist*MGeomCam.fConvMm2Deg";
|
---|
105 |
|
---|
106 | // --------------------------------------------------------------------------
|
---|
107 | //
|
---|
108 | // Display a result histogram either vs. size or energy
|
---|
109 | // FIXME: Could be moved into a new histogram class.
|
---|
110 | //
|
---|
111 | void MJTrainDisp::DisplayHist(TCanvas &c, Int_t i, MH3 &mh3) const
|
---|
112 | {
|
---|
113 | MH::SetPalette("pretty");
|
---|
114 |
|
---|
115 | TH1 &hist = *(TH1*)mh3.GetHist().Clone();
|
---|
116 | hist.SetBit(TH1::kNoStats);
|
---|
117 | hist.SetDirectory(0);
|
---|
118 |
|
---|
119 | TLine line;
|
---|
120 | line.SetLineStyle(kDashed);
|
---|
121 | line.SetLineWidth(1);
|
---|
122 |
|
---|
123 | c.cd(i+4);
|
---|
124 | gPad->SetBorderMode(0);
|
---|
125 | gPad->SetFrameBorderMode(0);
|
---|
126 | //gPad->SetFillColor(kWhite);
|
---|
127 | gPad->SetLogx();
|
---|
128 | gPad->SetGridx();
|
---|
129 | gPad->SetGridy();
|
---|
130 | //gPad->SetLeftMargin(0.12);
|
---|
131 | //gPad->SetRightMargin(0.12);
|
---|
132 |
|
---|
133 | const Float_t cutval = hist.GetYaxis()->GetBinLowEdge(4);
|
---|
134 |
|
---|
135 | TH1D heff;
|
---|
136 | heff.SetName(Form("Eff%s", hist.GetName()));
|
---|
137 | heff.SetTitle(Form("Cut efficiency vs. %s for \\vartheta<%.3f", hist.GetName(), TMath::Sqrt(cutval)));
|
---|
138 | heff.SetDirectory(0);
|
---|
139 | heff.SetXTitle(hist.GetXaxis()->GetTitle());
|
---|
140 | heff.SetYTitle("Efficiency");
|
---|
141 |
|
---|
142 | MH::SetBinning(&heff, hist.GetXaxis());
|
---|
143 |
|
---|
144 |
|
---|
145 | for (int x=0; x<=hist.GetNbinsX()+1; x++)
|
---|
146 | {
|
---|
147 | const Double_t n0 = hist.Integral(x, x, -1, -1);
|
---|
148 | if (n0>0)
|
---|
149 | heff.SetBinContent(x, hist.Integral(x, x, -1, 3)/n0);
|
---|
150 | }
|
---|
151 |
|
---|
152 | heff.SetMinimum(0);
|
---|
153 | heff.SetMaximum(1);
|
---|
154 | heff.DrawCopy();
|
---|
155 |
|
---|
156 | line.DrawLine(10, 0.5, 31623, 0.5);
|
---|
157 |
|
---|
158 | c.cd(i+0);
|
---|
159 | gPad->SetBorderMode(0);
|
---|
160 | gPad->SetFrameBorderMode(0);
|
---|
161 | //gPad->SetFillColor(kWhite);
|
---|
162 | gPad->SetLogx();
|
---|
163 | gPad->SetGridx();
|
---|
164 | gPad->SetGridy();
|
---|
165 | //gPad->SetLeftMargin(0.12);
|
---|
166 | //gPad->SetRightMargin(0.12);
|
---|
167 |
|
---|
168 | for (int x=0; x<=hist.GetNbinsX(); x++)
|
---|
169 | {
|
---|
170 | Float_t n = 0;
|
---|
171 | for (int y=1; y<=2; y++)
|
---|
172 | n += hist.GetBinContent(x,y);
|
---|
173 |
|
---|
174 | if (n==0)
|
---|
175 | continue;
|
---|
176 |
|
---|
177 | for (int y=0; y<=hist.GetNbinsY(); y++)
|
---|
178 | hist.SetBinContent(x, y, 200*hist.GetBinContent(x,y)/n);
|
---|
179 | }
|
---|
180 |
|
---|
181 | hist.SetMaximum(100);
|
---|
182 | hist.DrawCopy("colz");
|
---|
183 |
|
---|
184 | line.DrawLine(10, 0.04, 31623, 0.04);
|
---|
185 |
|
---|
186 | c.cd(i+2);
|
---|
187 | gPad->SetBorderMode(0);
|
---|
188 | gPad->SetFrameBorderMode(0);
|
---|
189 | //gPad->SetFillColor(kWhite);
|
---|
190 | gPad->SetLogx();
|
---|
191 | gPad->SetGridx();
|
---|
192 | gPad->SetGridy();
|
---|
193 | //gPad->SetLeftMargin(0.12);
|
---|
194 | //gPad->SetRightMargin(0.12);
|
---|
195 |
|
---|
196 | for (int x=0; x<=hist.GetNbinsX(); x++)
|
---|
197 | {
|
---|
198 | Float_t n = 0;
|
---|
199 | for (int y=0; y<=hist.GetNbinsY(); y++)
|
---|
200 | n += hist.GetBinContent(x,y);
|
---|
201 |
|
---|
202 | if (n==0)
|
---|
203 | continue;
|
---|
204 |
|
---|
205 | for (int y=0; y<=hist.GetNbinsY(); y++)
|
---|
206 | hist.SetBinContent(x, y, 100*hist.GetBinContent(x,y)/n);
|
---|
207 | }
|
---|
208 |
|
---|
209 | hist.SetMaximum(25);
|
---|
210 | hist.DrawCopy("colz");
|
---|
211 |
|
---|
212 | line.DrawLine(10, 0.04, 31623, 0.04);
|
---|
213 | }
|
---|
214 |
|
---|
215 | // --------------------------------------------------------------------------
|
---|
216 | //
|
---|
217 | // Display the result histograms in a new tab.
|
---|
218 | //
|
---|
219 | void MJTrainDisp::DisplayResult(MH3 &hsize, MH3 &henergy)
|
---|
220 | {
|
---|
221 | TCanvas &c = fDisplay->AddTab("Disp");
|
---|
222 | c.Divide(2,3);
|
---|
223 |
|
---|
224 | DisplayHist(c, 1, hsize);
|
---|
225 | DisplayHist(c, 2, henergy);
|
---|
226 | }
|
---|
227 |
|
---|
228 | // --------------------------------------------------------------------------
|
---|
229 | //
|
---|
230 | // Run Disp optimization
|
---|
231 | //
|
---|
232 | Bool_t MJTrainDisp::Train(const char *out, const MDataSet &set, Int_t num)
|
---|
233 | {
|
---|
234 | SetTitle(Form("Train%s: %s", fNameOutput.Data(), out));
|
---|
235 |
|
---|
236 | if (fDisplay)
|
---|
237 | fDisplay->SetTitle(fTitle);
|
---|
238 |
|
---|
239 | if (!set.IsValid())
|
---|
240 | {
|
---|
241 | *fLog << err << "ERROR - DataSet invalid!" << endl;
|
---|
242 | return kFALSE;
|
---|
243 | }
|
---|
244 |
|
---|
245 | if (!HasWritePermission(out))
|
---|
246 | return kFALSE;
|
---|
247 |
|
---|
248 | *fLog << inf;
|
---|
249 | fLog->Separator(GetDescriptor());
|
---|
250 |
|
---|
251 | // --------------------- Setup files --------------------
|
---|
252 | MReadMarsFile readtrn("Events");
|
---|
253 | MReadMarsFile readtst("Events");
|
---|
254 | readtrn.DisableAutoScheme();
|
---|
255 | readtst.DisableAutoScheme();
|
---|
256 |
|
---|
257 | if (!set.AddFilesOn(readtrn))
|
---|
258 | {
|
---|
259 | *fLog << err << "ERROR - Adding SequencesOn." << endl;
|
---|
260 | return kFALSE;
|
---|
261 | }
|
---|
262 | if (!set.AddFilesOff(readtst))
|
---|
263 | {
|
---|
264 | *fLog << err << "ERROR - Adding SequencesOff." << endl;
|
---|
265 | return kFALSE;
|
---|
266 | }
|
---|
267 |
|
---|
268 | // ----------------------- Setup Matrix ------------------
|
---|
269 | MHMatrix train("Train");
|
---|
270 | train.AddColumns(fRules);
|
---|
271 | if (fEnableWeights)
|
---|
272 | train.AddColumn("MWeight.fVal");
|
---|
273 | train.AddColumn(fTrainParameter);
|
---|
274 |
|
---|
275 | // ----------------------- Fill Matrix RF ----------------------
|
---|
276 | MTFillMatrix fill(fTitle);
|
---|
277 | fill.SetDisplay(fDisplay);
|
---|
278 | fill.SetLogStream(fLog);
|
---|
279 | fill.SetDestMatrix1(&train, num);
|
---|
280 | fill.SetReader(&readtrn);
|
---|
281 | fill.AddPreCuts(fPreCuts);
|
---|
282 | fill.AddPreCuts(fTrainCuts);
|
---|
283 | fill.AddPreTasks(fPreTasks);
|
---|
284 | fill.AddPostTasks(fPostTasks);
|
---|
285 | if (!fill.Process())
|
---|
286 | return kFALSE;
|
---|
287 |
|
---|
288 | // ------------------------ Train RF --------------------------
|
---|
289 | MRanForestCalc rf("Train", fTitle);
|
---|
290 | rf.SetNumTrees(fNumTrees);
|
---|
291 | rf.SetNdSize(fNdSize);
|
---|
292 | rf.SetNumTry(fNumTry);
|
---|
293 | rf.SetNumObsoleteVariables(1);
|
---|
294 | rf.SetLastDataColumnHasWeights(fEnableWeights);
|
---|
295 | rf.SetDisplay(fDisplay);
|
---|
296 | rf.SetLogStream(fLog);
|
---|
297 | rf.SetFileName(out);
|
---|
298 | rf.SetDebug(fDebug>1);
|
---|
299 | rf.SetNameOutput(fNameOutput);
|
---|
300 | rf.SetFunction(fResultFunction);
|
---|
301 |
|
---|
302 | /*
|
---|
303 | MBinning b(32, 10, 100000, "BinningEnergyEst", "log");
|
---|
304 |
|
---|
305 | if (!rf.TrainMultiRF(train, b.GetEdgesD())) // classification with one tree per bin
|
---|
306 | return;
|
---|
307 |
|
---|
308 | if (!rf.TrainSingleRF(train, b.GetEdgesD())) // classification into different bins
|
---|
309 | return;
|
---|
310 | */
|
---|
311 | if (!rf.TrainRegression(train)) // regression (best choice)
|
---|
312 | return kFALSE;
|
---|
313 |
|
---|
314 | // --------------------- Display result ----------------------
|
---|
315 |
|
---|
316 | gLog.Separator("Test");
|
---|
317 |
|
---|
318 | MH::SetPalette("pretty");
|
---|
319 |
|
---|
320 | MParList plist;
|
---|
321 | MTaskList tlist;
|
---|
322 | plist.AddToList(this);
|
---|
323 | plist.AddToList(&tlist);
|
---|
324 | //plist.AddToList(&b);
|
---|
325 |
|
---|
326 | MParameterD par("ThetaSquaredCut");
|
---|
327 | par.SetVal(fThetaCut*fThetaCut);
|
---|
328 | plist.AddToList(&par);
|
---|
329 |
|
---|
330 | MAlphaFitter fit;
|
---|
331 | fit.SetPolynomOrder(0);
|
---|
332 | fit.SetSignalFitMax(0.8);
|
---|
333 | fit.EnableBackgroundFit(kFALSE);
|
---|
334 | fit.SetSignalFunction(MAlphaFitter::kThetaSq);
|
---|
335 | fit.SetMinimizationStrategy(MAlphaFitter::kGaussSigma);
|
---|
336 | plist.AddToList(&fit);
|
---|
337 |
|
---|
338 | MFilterList list;
|
---|
339 | if (!list.AddToList(fPreCuts))
|
---|
340 | *fLog << err << "ERROR - Calling MFilterList::AddToList for fPreCuts failed!" << endl;
|
---|
341 | if (!list.AddToList(fTestCuts))
|
---|
342 | *fLog << err << "ERROR - Calling MFilterList::AddToList for fTestCuts failed!" << endl;
|
---|
343 |
|
---|
344 | MContinue cont(&list);
|
---|
345 | cont.SetInverted();
|
---|
346 |
|
---|
347 | const char *rule = "(MHillasSrc.fDist*MGeomCam.fConvMm2Deg)^2 + (Disp.fVal)^2 - (2*MHillasSrc.fDist*MGeomCam.fConvMm2Deg*Disp.fVal*cos(MHillasSrc.fAlpha*kDeg2Rad))";
|
---|
348 |
|
---|
349 | MParameterCalc calcthetasq(rule, "MThetaSqCalc");
|
---|
350 | calcthetasq.SetNameParameter("ThetaSquared");
|
---|
351 |
|
---|
352 | MChisqEval eval;
|
---|
353 | eval.SetY1("sqrt(ThetaSquared.fVal)");
|
---|
354 |
|
---|
355 | // ----------- Setup binnings ----------------
|
---|
356 | MBinning binsS(50, 10, 100000, "BinningSize", "log");
|
---|
357 | MBinning binsE(70, 10, 31623, "BinningEnergy", "log");
|
---|
358 | MBinning binsG(50, -10, 10, "BinningSlope", "lin");
|
---|
359 | MBinning binsX(50, -1, 1, "BinningResidualDist", "lin");
|
---|
360 | MBinning binsL(50, 0, 0.3, "BinningLeakage", "lin");
|
---|
361 | MBinning binsT(51, -0.005, 0.505, "BinningTheta", "asin");
|
---|
362 | MBinning binsC(50, 1e-2, 1, "BinningConc", "log");
|
---|
363 | MBinning binsW(50, 0, 0.5, "BinningLength", "lin");
|
---|
364 | MBinning binsM(50, 0, 0.3, "BinningWidth", "lin");
|
---|
365 | MBinning binsV(75, 0, par.GetVal()*25, "BinningThetaSq", "lin");
|
---|
366 |
|
---|
367 | plist.AddToList(&binsG);
|
---|
368 | plist.AddToList(&binsS);
|
---|
369 | plist.AddToList(&binsX);
|
---|
370 | plist.AddToList(&binsE);
|
---|
371 | plist.AddToList(&binsL);
|
---|
372 | plist.AddToList(&binsT);
|
---|
373 | plist.AddToList(&binsC);
|
---|
374 | plist.AddToList(&binsV);
|
---|
375 | plist.AddToList(&binsW);
|
---|
376 | plist.AddToList(&binsM);
|
---|
377 |
|
---|
378 | // ----------- Setup some histograms ----------------
|
---|
379 |
|
---|
380 | MHThetaSq hist;
|
---|
381 | hist.SkipHistTime();
|
---|
382 | hist.SkipHistTheta();
|
---|
383 | hist.SkipHistEnergy();
|
---|
384 |
|
---|
385 | // To speed it up we could precalculate it.
|
---|
386 | const char *res = "Disp.fVal-MHillasSrc.fDist*3.37e-3";
|
---|
387 |
|
---|
388 | MHn hres1("Disp1", "Xi Residual (Dist/Disp)");
|
---|
389 | hres1.AddHist("MHillas.fSize", res);
|
---|
390 | hres1.InitName("ResSize;Size;ResidualDist");
|
---|
391 | hres1.InitTitle(";S [phe];Disp-Dist [\\circ];");
|
---|
392 | hres1.SetDrawOption("colz profx");
|
---|
393 | hres1.AddHist("MHillasExt.fSlopeLong*sign(MHillasSrc.fCosDeltaAlpha)/3.37e-3", res);
|
---|
394 | hres1.InitName("ResSlope;Slope;ResidualDist");
|
---|
395 | hres1.InitTitle(";Slope;Disp-Dist [\\circ];");
|
---|
396 | hres1.SetDrawOption("colz profx");
|
---|
397 | hres1.AddHist("MNewImagePar.fLeakage1", res);
|
---|
398 | hres1.InitName("ResLeak;Leakage;ResidualDist");
|
---|
399 | hres1.InitTitle(";Leak;Disp-Dist [\\circ];");
|
---|
400 | hres1.SetDrawOption("colz profx");
|
---|
401 | hres1.AddHist("MPointingPos.fZd", res);
|
---|
402 | hres1.InitName("ResTheta;Theta;ResidualDist");
|
---|
403 | hres1.InitTitle(";Zd [\\circ];Disp-Dist [\\circ];");
|
---|
404 | hres1.SetDrawOption("colz profx");
|
---|
405 |
|
---|
406 | MHn hres2("Disp2", "Dist Residual (Disp-Dist)");
|
---|
407 | hres2.AddHist("MHillas.fLength*3.37e-3", res);
|
---|
408 | hres2.InitName("ResLength;Length;ResidualDist");
|
---|
409 | hres2.InitTitle(";L [\\circ];Disp-Dist [\\circ];");
|
---|
410 | hres2.SetDrawOption("colz profx");
|
---|
411 | hres2.AddHist("MNewImagePar.fConc1", res);
|
---|
412 | hres2.InitName("ResConc1;Conc;ResidualDist");
|
---|
413 | hres2.InitTitle(";C;Disp-Dist [\\circ];");
|
---|
414 | hres2.SetDrawOption("colz profx");
|
---|
415 | hres2.AddHist("MHillas.fWidth*3.37e-3", res);
|
---|
416 | hres2.InitName("ResWidth;Width;ResidualDist");
|
---|
417 | hres2.InitTitle(";W [\\circ];Disp-Dist [\\circ];");
|
---|
418 | hres2.SetDrawOption("colz profx");
|
---|
419 | hres2.AddHist("MMcEvt.fEnergy", res);
|
---|
420 | hres2.InitName("ResEmc;Energy;ResidualDist");
|
---|
421 | hres2.InitTitle(";E_{mc} [GeV];Disp-Dist [\\circ];");
|
---|
422 | hres2.SetDrawOption("colz profx");
|
---|
423 |
|
---|
424 | MH3 hdisp1("MHillas.fSize", "ThetaSquared.fVal");
|
---|
425 | MH3 hdisp2("MMcEvt.fEnergy", "ThetaSquared.fVal");
|
---|
426 | hdisp1.SetName("Size;Size;ThetaSq");
|
---|
427 | hdisp2.SetName("Energy;Energy;ThetaSq");
|
---|
428 | hdisp1.SetTitle("\\vartheta distribution vs. Size:Size [phe]:\\vartheta^2 [\\circ]");
|
---|
429 | hdisp2.SetTitle("\\vartheta distribution vs. Energy:Energy [GeV]:\\vartheta^2 [\\circ]");
|
---|
430 |
|
---|
431 | // -------------- Setup fill tasks ----------------
|
---|
432 |
|
---|
433 | MFillH fillh(&hist, "", "FillThetaSq");
|
---|
434 | MFillH fillh2a(&hres1, "", "FillResiduals1");
|
---|
435 | MFillH fillh2b(&hres2, "", "FillResiduals2");
|
---|
436 | MFillH fillh2c(&hdisp1, "", "FillSize");
|
---|
437 | MFillH fillh2d(&hdisp2, "", "FillEnergy");
|
---|
438 | fillh2c.SetBit(MFillH::kDoNotDisplay);
|
---|
439 | fillh2d.SetBit(MFillH::kDoNotDisplay);
|
---|
440 |
|
---|
441 | // --------------- Setup weighting -------------------
|
---|
442 |
|
---|
443 | if (fEnableWeights)
|
---|
444 | {
|
---|
445 | fillh.SetWeight();
|
---|
446 | fillh2a.SetWeight();
|
---|
447 | fillh2b.SetWeight();
|
---|
448 | fillh2c.SetWeight();
|
---|
449 | fillh2d.SetWeight();
|
---|
450 | eval.SetNameWeight();
|
---|
451 | }
|
---|
452 |
|
---|
453 | // --------------- Setup tasklist -------------------
|
---|
454 |
|
---|
455 | tlist.AddToList(&readtst);
|
---|
456 | tlist.AddToList(fPreTasks);
|
---|
457 | tlist.AddToList(&cont);
|
---|
458 | tlist.AddToList(&rf);
|
---|
459 | tlist.AddToList(&calcthetasq);
|
---|
460 | tlist.AddToList(fPostTasks);
|
---|
461 | tlist.AddToList(&fillh);
|
---|
462 | tlist.AddToList(&fillh2a);
|
---|
463 | tlist.AddToList(&fillh2b);
|
---|
464 | tlist.AddToList(&fillh2c);
|
---|
465 | tlist.AddToList(&fillh2d);
|
---|
466 | tlist.AddToList(fTestTasks);
|
---|
467 | tlist.AddToList(&eval);
|
---|
468 |
|
---|
469 | // ------------- Setup/run eventloop -----------------
|
---|
470 |
|
---|
471 | MEvtLoop loop(fTitle);
|
---|
472 | loop.SetLogStream(fLog);
|
---|
473 | loop.SetDisplay(fDisplay);
|
---|
474 | loop.SetParList(&plist);
|
---|
475 | //if (!SetupEnv(loop))
|
---|
476 | // return kFALSE;
|
---|
477 |
|
---|
478 | if (!loop.Eventloop())
|
---|
479 | return kFALSE;
|
---|
480 |
|
---|
481 | // ---------------- Prepare result -------------------
|
---|
482 |
|
---|
483 | // Print the result
|
---|
484 | *fLog << inf;
|
---|
485 | *fLog << "Rule: " << rule << endl;
|
---|
486 | *fLog << "Disp: " << fTrainParameter << endl;
|
---|
487 | hist.GetAlphaFitter().Print("result");
|
---|
488 |
|
---|
489 | // The user has closed the display
|
---|
490 | if (!fDisplay)
|
---|
491 | return kTRUE;
|
---|
492 |
|
---|
493 | DisplayResult(hdisp1, hdisp2);
|
---|
494 |
|
---|
495 | SetPathOut(out);
|
---|
496 | return WriteDisplay(0, "UPDATE");
|
---|
497 | }
|
---|
498 |
|
---|
499 | /*
|
---|
500 | #include "MParameterCalc.h"
|
---|
501 | #include "MHillasCalc.h"
|
---|
502 | #include "../mpointing/MSrcPosRndm.h"
|
---|
503 |
|
---|
504 | Bool_t MJTrainDisp::TrainGhostbuster(const char *out, const MDataSet &set, Int_t num)
|
---|
505 | {
|
---|
506 | SetTitle(Form("TrainGhostbuster: %s", out));
|
---|
507 |
|
---|
508 | if (fDisplay)
|
---|
509 | fDisplay->SetTitle(fTitle);
|
---|
510 |
|
---|
511 | if (!set.IsValid())
|
---|
512 | {
|
---|
513 | *fLog << err << "ERROR - DataSet invalid!" << endl;
|
---|
514 | return kFALSE;
|
---|
515 | }
|
---|
516 |
|
---|
517 | if (!HasWritePermission(out))
|
---|
518 | return kFALSE;
|
---|
519 |
|
---|
520 | *fLog << inf;
|
---|
521 | fLog->Separator(GetDescriptor());
|
---|
522 |
|
---|
523 | // --------------------- Setup files --------------------
|
---|
524 | MReadMarsFile readtrn("Events");
|
---|
525 | MReadMarsFile readtst("Events");
|
---|
526 | readtrn.DisableAutoScheme();
|
---|
527 | readtst.DisableAutoScheme();
|
---|
528 |
|
---|
529 | if (!set.AddFilesOn(readtrn))
|
---|
530 | return kFALSE;
|
---|
531 | if (!set.AddFilesOff(readtst))
|
---|
532 | return kFALSE;
|
---|
533 |
|
---|
534 | // ----------------------- Setup Matrix ------------------
|
---|
535 | MHMatrix train("Train");
|
---|
536 | train.AddColumns(fRules);
|
---|
537 | if (fEnableWeights)
|
---|
538 | train.AddColumn("MWeight.fVal");
|
---|
539 | train.AddColumn("sign(MHillasSrc.fCosDeltaAlpha)==sign(SignStore.fVal)");
|
---|
540 |
|
---|
541 | MParameterCalc calc("MHillasSrc.fCosDeltaAlpha", "SignStore");
|
---|
542 | calc.SetNameParameter("SignStore");
|
---|
543 |
|
---|
544 | MSrcPosRndm rndm;
|
---|
545 | rndm.SetRule(fTrainParameter);
|
---|
546 | //rndm.SetDistOfSource(120*3.37e-3);
|
---|
547 |
|
---|
548 | MHillasCalc hcalc;
|
---|
549 | hcalc.SetFlags(MHillasCalc::kCalcHillasSrc);
|
---|
550 |
|
---|
551 | // ----------------------- Fill Matrix RF ----------------------
|
---|
552 | MTFillMatrix fill(fTitle);
|
---|
553 | fill.SetDisplay(fDisplay);
|
---|
554 | fill.SetLogStream(fLog);
|
---|
555 | fill.SetDestMatrix1(&train, num);
|
---|
556 | fill.SetReader(&readtrn);
|
---|
557 | fill.AddPreCuts(fPreCuts);
|
---|
558 | fill.AddPreCuts(fTrainCuts);
|
---|
559 | fill.AddPreTasks(fPreTasks);
|
---|
560 | fill.AddPostTasks(fPostTasks);
|
---|
561 | fill.AddPostTask(&calc);
|
---|
562 | fill.AddPostTask(&rndm);
|
---|
563 | fill.AddPostTask(&hcalc);
|
---|
564 | if (!fill.Process())
|
---|
565 | return kFALSE;
|
---|
566 |
|
---|
567 | // ------------------------ Train RF --------------------------
|
---|
568 | MRanForestCalc rf("TrainGhostbuster", fTitle);
|
---|
569 | rf.SetNumTrees(fNumTrees);
|
---|
570 | rf.SetNdSize(fNdSize);
|
---|
571 | rf.SetNumTry(fNumTry);
|
---|
572 | rf.SetNumObsoleteVariables(1);
|
---|
573 | rf.SetLastDataColumnHasWeights(fEnableWeights);
|
---|
574 | rf.SetDisplay(fDisplay);
|
---|
575 | rf.SetLogStream(fLog);
|
---|
576 | rf.SetFileName(out);
|
---|
577 | rf.SetDebug(fDebug>1);
|
---|
578 | rf.SetNameOutput("Sign");
|
---|
579 |
|
---|
580 | if (!rf.TrainRegression(train)) // regression (best choice)
|
---|
581 | return kFALSE;
|
---|
582 |
|
---|
583 | // --------------------- Display result ----------------------
|
---|
584 |
|
---|
585 | gLog.Separator("Test");
|
---|
586 |
|
---|
587 | MH::SetPalette("pretty");
|
---|
588 |
|
---|
589 | MParList plist;
|
---|
590 | MTaskList tlist;
|
---|
591 | plist.AddToList(this);
|
---|
592 | plist.AddToList(&tlist);
|
---|
593 | //plist.AddToList(&b);
|
---|
594 |
|
---|
595 | MFilterList list;
|
---|
596 | if (!list.AddToList(fPreCuts))
|
---|
597 | *fLog << err << "ERROR - Calling MFilterList::AddToList for fPreCuts failed!" << endl;
|
---|
598 | if (!list.AddToList(fTestCuts))
|
---|
599 | *fLog << err << "ERROR - Calling MFilterList::AddToList for fTestCuts failed!" << endl;
|
---|
600 |
|
---|
601 | MContinue cont(&list);
|
---|
602 | cont.SetInverted();
|
---|
603 |
|
---|
604 | const char *rule = "abs(Sign.fVal-(sign(MHillasSrc.fCosDeltaAlpha)==sign(SignStore.fVal)))";
|
---|
605 |
|
---|
606 | //MChisqEval eval;
|
---|
607 | //eval.SetY1("sqrt(ThetaSquared.fVal)");
|
---|
608 |
|
---|
609 | MH3 hsign1("MHillas.fSize", rule);
|
---|
610 | MH3 hsign2("MMcEvt.fEnergy", rule);
|
---|
611 | hsign1.SetTitle("Was ist das? vs. Size:Size [phe]:XXX");
|
---|
612 | hsign2.SetTitle("Was ist das? vs. Energy:Enerhy [GeV]:XXXX");
|
---|
613 |
|
---|
614 | MBinning binsx( 70, 10, 31623, "BinningMH3X", "log");
|
---|
615 | MBinning binsy( 51, 0, 1, "BinningMH3Y", "lin");
|
---|
616 |
|
---|
617 | plist.AddToList(&binsx);
|
---|
618 | plist.AddToList(&binsy);
|
---|
619 |
|
---|
620 | MFillH fillh2a(&hsign1, "", "FillSize");
|
---|
621 | MFillH fillh2b(&hsign2, "", "FillEnergy");
|
---|
622 | fillh2a.SetDrawOption("profx colz");
|
---|
623 | fillh2b.SetDrawOption("profx colz");
|
---|
624 | fillh2a.SetNameTab("Size");
|
---|
625 | fillh2b.SetNameTab("Energy");
|
---|
626 |
|
---|
627 | tlist.AddToList(&readtst);
|
---|
628 | tlist.AddToList(&cont);
|
---|
629 | tlist.AddToList(&calc);
|
---|
630 | tlist.AddToList(&rndm);
|
---|
631 | tlist.AddToList(&hcalc);
|
---|
632 | tlist.AddToList(&rf);
|
---|
633 | tlist.AddToList(&fillh2a);
|
---|
634 | tlist.AddToList(&fillh2b);
|
---|
635 | //tlist.AddToList(&eval);
|
---|
636 |
|
---|
637 | MEvtLoop loop(fTitle);
|
---|
638 | loop.SetLogStream(fLog);
|
---|
639 | loop.SetDisplay(fDisplay);
|
---|
640 | loop.SetParList(&plist);
|
---|
641 | //if (!SetupEnv(loop))
|
---|
642 | // return kFALSE;
|
---|
643 |
|
---|
644 | if (!loop.Eventloop())
|
---|
645 | return kFALSE;
|
---|
646 |
|
---|
647 | // Print the result
|
---|
648 | *fLog << inf << "Rule: " << rule << endl;
|
---|
649 |
|
---|
650 | //DisplayResult(hdisp1, hdisp2);
|
---|
651 |
|
---|
652 | SetPathOut(out);
|
---|
653 | if (!WriteDisplay(0, "UPDATE"))
|
---|
654 | return kFALSE;
|
---|
655 |
|
---|
656 | return kTRUE;
|
---|
657 | }
|
---|
658 | */
|
---|