source: branches/AddingGoogleTestEnvironment/mjtrain/MJTrainImpact.cc@ 19779

Last change on this file since 19779 was 9846, checked in by tbretz, 14 years ago
Added new class MJTrainImpact to estimate the impact parameter by use of a random forest.
File size: 10.5 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 8/2010 <mailto:thomas.bretz@epfl.ch>
19!
20! Copyright: MAGIC Software Development, 2010
21!
22!
23\* ======================================================================== */
24
25/////////////////////////////////////////////////////////////////////////////
26//
27// MJTrainImpact
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// MJTrainImpact opt;
37// //opt.SetDebug();
38// opt.AddParameter("MHillas.fSize");
39// opt.AddParameter("MHillasSrc.fDist");
40// MStatusDisplay *d = new MStatusDisplay;
41// opt.SetDisplay(d);
42// opt.AddPreCut("MHillasSrc.fDCA*MGeomCam.fConvMm2Deg<0.3");
43// opt.Train("rf-impact.root", set, 30000); // Number of train events
44//
45// Random Numbers:
46// ---------------
47// Use:
48// if(gRandom)
49// delete gRandom;
50// gRandom = new TRandom3();
51// in advance to change the random number generator.
52//
53////////////////////////////////////////////////////////////////////////////
54#include "MJTrainImpact.h"
55
56#include "MHMatrix.h"
57
58#include "MLog.h"
59#include "MLogManip.h"
60
61// tools
62#include "MDataSet.h"
63#include "MTFillMatrix.h"
64#include "MStatusDisplay.h"
65
66// eventloop
67#include "MParList.h"
68#include "MTaskList.h"
69#include "MEvtLoop.h"
70
71// tasks
72#include "MReadMarsFile.h"
73#include "MContinue.h"
74#include "MFillH.h"
75#include "MRanForestCalc.h"
76
77// histograms
78#include "MHn.h"
79#include "MBinning.h"
80#include "MHEnergyEst.h"
81
82// filter
83#include "MFilterList.h"
84
85ClassImp(MJTrainImpact);
86
87using namespace std;
88
89Bool_t MJTrainImpact::Train(const char *out, const MDataSet &set, Int_t num)
90{
91 SetTitle(Form("Train%s: %s", fNameOutput.Data(), out));
92
93 if (fDisplay)
94 fDisplay->SetTitle(fTitle);
95
96 if (!set.IsValid())
97 {
98 *fLog << err << "ERROR - DataSet invalid!" << endl;
99 return kFALSE;
100 }
101
102 if (!HasWritePermission(out))
103 return kFALSE;
104
105 *fLog << inf;
106 fLog->Separator(GetDescriptor());
107
108 // --------------------- Setup files --------------------
109 MReadMarsFile readtrn("Events");
110 MReadMarsFile readtst("Events");
111 readtrn.DisableAutoScheme();
112 readtst.DisableAutoScheme();
113
114 if (!set.AddFilesOn(readtrn))
115 {
116 *fLog << err << "ERROR - Adding SequencesOn." << endl;
117 return kFALSE;
118 }
119 if (!set.AddFilesOff(readtst))
120 {
121 *fLog << err << "ERROR - Adding SequencesOff." << endl;
122 return kFALSE;
123 }
124
125 // ----------------------- Setup Matrix ------------------
126 MHMatrix train("Train");
127 train.AddColumns(fRules);
128 if (fEnableWeights)
129 train.AddColumn("MWeight.fVal");
130 train.AddColumn(fTrainParameter);
131
132 // ----------------------- Fill Matrix RF ----------------------
133 MTFillMatrix fill(fTitle);
134 fill.SetDisplay(fDisplay);
135 fill.SetLogStream(fLog);
136 fill.SetDestMatrix1(&train, num);
137 fill.SetReader(&readtrn);
138 fill.AddPreCuts(fPreCuts);
139 fill.AddPreCuts(fTrainCuts);
140 fill.AddPreTasks(fPreTasks);
141 fill.AddPostTasks(fPostTasks);
142 if (!fill.Process())
143 return kFALSE;
144
145 // ------------------------ Train RF --------------------------
146 MRanForestCalc rf("Train", fTitle);
147 rf.SetNumTrees(fNumTrees);
148 rf.SetNdSize(fNdSize);
149 rf.SetNumTry(fNumTry);
150 rf.SetNumObsoleteVariables(1);
151 rf.SetLastDataColumnHasWeights(fEnableWeights);
152 rf.SetDisplay(fDisplay);
153 rf.SetLogStream(fLog);
154 rf.SetFileName(out);
155 rf.SetDebug(fDebug>1);
156 rf.SetNameOutput(fNameOutput);
157 rf.SetFunction(fResultFunction);
158
159 /*
160 MBinning b(32, 10, 100000, "BinningEnergyEst", "log");
161
162 if (!rf.TrainMultiRF(train, b.GetEdgesD())) // classification with one tree per bin
163 return;
164
165 if (!rf.TrainSingleRF(train, b.GetEdgesD())) // classification into different bins
166 return;
167 */
168 if (!rf.TrainRegression(train)) // regression (best choice)
169 return kFALSE;
170
171 // --------------------- Display result ----------------------
172
173 gLog.Separator("Test");
174
175 MH::SetPalette("pretty");
176
177 MParList plist;
178 MTaskList tlist;
179 plist.AddToList(this);
180 plist.AddToList(&tlist);
181 //plist.AddToList(&b);
182
183 MFilterList list;
184 if (!list.AddToList(fPreCuts))
185 *fLog << err << "ERROR - Calling MFilterList::AddToList for fPreCuts failed!" << endl;
186 if (!list.AddToList(fTestCuts))
187 *fLog << err << "ERROR - Calling MFilterList::AddToList for fTestCuts failed!" << endl;
188
189 MContinue cont(&list);
190 cont.SetInverted();
191
192 // -------------------------------------------------------------
193 MBinning binsS(40, 10, 100000, "BinningSize", "log");
194 MBinning binsE(40, 10, 100000, "BinningEnergy", "log");
195// MBinning binsF(40, 10, 100000, "BinningEnergyEst", "log");
196 MBinning binsG(90, -14, 14, "BinningSlope", "lin");
197 MBinning binsR(50,-100, 100, "BinningImpactResidual", "lin");
198 MBinning binsL(50, 0, 0.3, "BinningLeakage", "lin");
199 MBinning binsT(51, -0.005, 0.505, "BinningTheta", "asin");
200 MBinning binsD(70, 0, 2.4, "BinningDist", "lin");
201 MBinning binsC(50, 1e-2, 1, "BinningConc", "log");
202 MBinning binsI(32, 0, 800, "BinningImpact", "lin");
203
204 plist.AddToList(&binsG);
205 plist.AddToList(&binsS);
206 plist.AddToList(&binsR);
207 plist.AddToList(&binsE);
208// plist.AddToList(&binsF);
209 plist.AddToList(&binsL);
210 plist.AddToList(&binsT);
211 plist.AddToList(&binsD);
212 plist.AddToList(&binsC);
213 plist.AddToList(&binsI);
214
215 //MHEnergyEst hist;
216
217 // To speed it up we could precalculate it.
218 const char *res = "(MImpactEst.fVal-MMcEvt.fImpact)/100";
219
220 MHn hres1("Impact1", "Impact Residual (I_{est} - I_{mc})");
221 hres1.AddHist("MHillas.fSize", res);
222 hres1.InitName("ResSize;Size;ImpactResidual");
223 hres1.InitTitle(";S [phe];\\Delta I;");
224 hres1.SetDrawOption("colz profx");
225 hres1.AddHist("MHillasExt.fSlopeLong*sign(MHillasSrc.fCosDeltaAlpha)/MGeomCam.fConvMm2Deg", res);
226 hres1.InitName("ResSlope;Slope;ImpactResidual");
227 hres1.InitTitle(";Slope;\\Delta I;");
228 hres1.SetDrawOption("colz profx");
229 hres1.AddHist("MNewImagePar.fLeakage1", res);
230 hres1.InitName("ResLeak;Leakage;ImpactResidual");
231 hres1.InitTitle(";Leak;\\Delta I;");
232 hres1.SetDrawOption("colz profx");
233 hres1.AddHist("MHillasSrc.fDist*MGeomCam.fConvMm2Deg", res);
234 hres1.InitName("ResDist;Dist;ImpactResidual");
235 hres1.InitTitle(";D [\\circ];\\Delta I;");
236 hres1.SetDrawOption("colz profx");
237
238 MHn hres2("Impact2", "Impact Residual (I_{est} - I_{mc})");
239 hres2.AddHist("MMcEvt.fEnergy", res);
240 hres2.InitName("ResEmc;Energy;ImpactResidual");
241 hres2.InitTitle(";E_{mc} [m];\\Delta I;");
242 hres2.SetDrawOption("colz profx");
243 hres2.SetAutoRange(kFALSE, kFALSE, kFALSE);
244 hres2.AddHist("MPointingPos.fZd", res);
245 hres2.InitName("ResTheta;Theta;ImpactResidual");
246 hres2.InitTitle(";Zd [\\circ];\\Delta I;");
247 hres2.SetDrawOption("colz profx");
248 hres2.AddHist("MImpactEst.fVal/100", res);
249 hres2.InitName("ResIest;Impact;ImpactResidual");
250 hres2.InitTitle(";I_{est} [m];\\Delta I;");
251 hres2.SetDrawOption("colz profx");
252 hres2.SetAutoRange(kFALSE, kFALSE, kFALSE);
253 hres2.AddHist("MNewImagePar.fConc1", res);
254 hres2.InitName("ResConc1;Conc;ImpactResidual");
255 hres2.InitTitle(";C;\\Delta I;");
256 hres2.SetDrawOption("colz profx");
257
258 MHn hres3("Resolution", "Impact Resolution");
259 hres3.AddHist("MImpactEst.fVal/100", "(MMcEvt.fImpact/MImpactEst.fVal-1)^2", MH3::kProfile);
260 hres3.InitName("ResIest;Impact;");
261 hres3.InitTitle(";I_{est} [m];Resolution \\sqrt{<(I_{mc}/I_{est}-1)^{2}>};");
262 hres3.SetConversion("sqrt(x)");
263
264 hres3.AddHist("MHillas.fSize", "(MMcEvt.fImpact/MImpactEst.fVal-1)^2", MH3::kProfile);
265 hres3.InitName("ResSize;Size;");
266 hres3.InitTitle(";S [phe];Resolution \\sqrt{<(I_{mc}/I_{est}-1)^{2}>};");
267 hres3.SetConversion("sqrt(x)");
268/*
269 hres3.AddHist("MMcEvt.fEnergy", "(MImpactEst.fVal/MMcEvt.fImpact-1)^2", MH3::kProfile);
270 hres3.InitName("ResEmc;EnergyEst;");
271 hres3.InitTitle(";E_{mc} [GeV];Resolution \\sqrt{<(I_{mc}/I_{est}-1)^{2}>};");
272 hres3.SetConversion("sqrt(x)");
273 */
274 hres3.AddHist("MMcEvt.fEnergy", "(MMcEvt.fImpact/MImpactEst.fVal-1)^2", MH3::kProfile);
275 hres3.InitName("ResEnergy;Energy;");
276 hres3.InitTitle(";E_{mc} [m];Resolution \\sqrt{<(I_{mc}/I_{est}-1)^{2}>};");
277 hres3.SetConversion("sqrt(x)");
278
279 hres3.AddHist("MPointingPos.fZd", "(MMcEvt.fImpact/MImpactEst.fVal-1)^2", MH3::kProfile);
280 hres3.InitName("ResTheta;Theta;");
281 hres3.InitTitle(";\\Theta [\\circ];Resolution \\sqrt{<(I_{mc}/I_{est}-1)^{2}>};");
282 hres3.SetConversion("sqrt(x)");
283
284 //MFillH fillh0(&hist);
285 MFillH fillh1(&hres1, "", "FillResiduals1");
286 MFillH fillh2(&hres2, "", "FillResiduals2");
287 MFillH fillh3(&hres3, "", "FillResolution");
288
289 if (fEnableWeights)
290 {
291 //fillh0.SetWeight();
292 fillh1.SetWeight();
293 fillh2.SetWeight();
294 fillh3.SetWeight();
295 }
296
297 tlist.AddToList(&readtst);
298 tlist.AddToList(fPreTasks);
299 tlist.AddToList(&cont);
300 tlist.AddToList(&rf);
301 tlist.AddToList(fPostTasks);
302 //tlist.AddToList(&fillh0);
303 tlist.AddToList(&fillh1);
304 tlist.AddToList(&fillh2);
305 tlist.AddToList(&fillh3);
306 tlist.AddToList(fTestTasks);
307
308 MEvtLoop loop(fTitle);
309 loop.SetLogStream(fLog);
310 loop.SetDisplay(fDisplay);
311 loop.SetParList(&plist);
312 //if (!SetupEnv(loop))
313 // return kFALSE;
314
315 if (!loop.Eventloop())
316 return kFALSE;
317
318 TObjArray arr;
319 arr.Add(const_cast<MDataSet*>(&set));
320 if (fDisplay)
321 arr.Add(fDisplay);
322
323 SetPathOut(out);
324 return WriteContainer(arr, 0, "UPDATE");
325}
Note: See TracBrowser for help on using the repository browser.