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