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, 2006
|
---|
21 | !
|
---|
22 | !
|
---|
23 | \* ======================================================================== */
|
---|
24 |
|
---|
25 | /////////////////////////////////////////////////////////////////////////////
|
---|
26 | //
|
---|
27 | // MJTrainSeparation
|
---|
28 | //
|
---|
29 | ////////////////////////////////////////////////////////////////////////////
|
---|
30 | #include "MJTrainSeparation.h"
|
---|
31 |
|
---|
32 | #include <TF1.h>
|
---|
33 | #include <TH2.h>
|
---|
34 | #include <TChain.h>
|
---|
35 | #include <TGraph.h>
|
---|
36 | #include <TCanvas.h>
|
---|
37 | #include <TVirtualPad.h>
|
---|
38 |
|
---|
39 | #include "MHMatrix.h"
|
---|
40 |
|
---|
41 | #include "MLog.h"
|
---|
42 | #include "MLogManip.h"
|
---|
43 |
|
---|
44 | // tools
|
---|
45 | #include "MMath.h"
|
---|
46 | #include "MDataSet.h"
|
---|
47 | #include "MTFillMatrix.h"
|
---|
48 | #include "MChisqEval.h"
|
---|
49 | #include "MStatusDisplay.h"
|
---|
50 |
|
---|
51 | // eventloop
|
---|
52 | #include "MParList.h"
|
---|
53 | #include "MTaskList.h"
|
---|
54 | #include "MEvtLoop.h"
|
---|
55 |
|
---|
56 | // tasks
|
---|
57 | #include "MReadMarsFile.h"
|
---|
58 | #include "MContinue.h"
|
---|
59 | #include "MFillH.h"
|
---|
60 | #include "MRanForestCalc.h"
|
---|
61 | #include "MParameterCalc.h"
|
---|
62 |
|
---|
63 | // container
|
---|
64 | #include "MMcEvt.hxx"
|
---|
65 | #include "MParameters.h"
|
---|
66 |
|
---|
67 | // histograms
|
---|
68 | #include "MBinning.h"
|
---|
69 | #include "MH3.h"
|
---|
70 | #include "MHHadronness.h"
|
---|
71 |
|
---|
72 | // filter
|
---|
73 | #include "MF.h"
|
---|
74 | #include "MFEventSelector.h"
|
---|
75 | #include "MFilterList.h"
|
---|
76 |
|
---|
77 | ClassImp(MJTrainSeparation);
|
---|
78 |
|
---|
79 | using namespace std;
|
---|
80 |
|
---|
81 | void MJTrainSeparation::DisplayResult(MH3 &h31, MH3 &h32)
|
---|
82 | {
|
---|
83 | TH2D &g = (TH2D&)h32.GetHist();
|
---|
84 | TH2D &h = (TH2D&)h31.GetHist();
|
---|
85 |
|
---|
86 | h.SetMarkerColor(kRed);
|
---|
87 | g.SetMarkerColor(kGreen);
|
---|
88 |
|
---|
89 | TH2D res1(g);
|
---|
90 | TH2D res2(g);
|
---|
91 |
|
---|
92 | h.SetTitle("Hadronness-Distribution vs. Size");
|
---|
93 | res1.SetTitle("Significance Li/Ma");
|
---|
94 | res1.SetXTitle("Size [phe]");
|
---|
95 | res1.SetYTitle("Hadronness");
|
---|
96 | res2.SetTitle("Significance-Distribution");
|
---|
97 | res2.SetXTitle("Size-Cut [phe]");
|
---|
98 | res2.SetYTitle("Hadronness-Cut");
|
---|
99 | res1.SetContour(50);
|
---|
100 | res2.SetContour(50);
|
---|
101 |
|
---|
102 | const Int_t nx = h.GetNbinsX();
|
---|
103 | const Int_t ny = h.GetNbinsY();
|
---|
104 |
|
---|
105 | gROOT->SetSelectedPad(NULL);
|
---|
106 |
|
---|
107 | TGraph gr1;
|
---|
108 | TGraph gr2;
|
---|
109 | for (int x=0; x<nx; x++)
|
---|
110 | {
|
---|
111 | TH1 *hx = h.ProjectionY("H_py", x+1, x+1);
|
---|
112 | TH1 *gx = g.ProjectionY("G_py", x+1, x+1);
|
---|
113 |
|
---|
114 | Double_t max1 = -1;
|
---|
115 | Double_t max2 = -1;
|
---|
116 | Int_t maxy1 = 0;
|
---|
117 | Int_t maxy2 = 0;
|
---|
118 | for (int y=0; y<ny; y++)
|
---|
119 | {
|
---|
120 | const Float_t s = gx->Integral(1, y+1);
|
---|
121 | const Float_t b = hx->Integral(1, y+1);
|
---|
122 | const Float_t sig1 = MMath::SignificanceLiMa(s+b, b);
|
---|
123 | const Float_t sig2 = s<1 ? 0 : MMath::SignificanceLiMa(s+b, b)*TMath::Log10(s);
|
---|
124 | if (sig1>max1)
|
---|
125 | {
|
---|
126 | maxy1 = y;
|
---|
127 | max1 = sig1;
|
---|
128 | }
|
---|
129 | if (sig2>max2)
|
---|
130 | {
|
---|
131 | maxy2 = y;
|
---|
132 | max2 = sig2;
|
---|
133 | }
|
---|
134 |
|
---|
135 | res1.SetBinContent(x+1, y+1, sig1);
|
---|
136 | }
|
---|
137 |
|
---|
138 | gr1.SetPoint(x, h.GetXaxis()->GetBinCenter(x+1), h.GetYaxis()->GetBinCenter(maxy1+1));
|
---|
139 | gr2.SetPoint(x, h.GetXaxis()->GetBinCenter(x+1), h.GetYaxis()->GetBinCenter(maxy2+1));
|
---|
140 |
|
---|
141 | delete hx;
|
---|
142 | delete gx;
|
---|
143 | }
|
---|
144 |
|
---|
145 | for (int x=0; x<nx; x++)
|
---|
146 | {
|
---|
147 | TH1 *hx = h.ProjectionY("H_py", x+1);
|
---|
148 | TH1 *gx = g.ProjectionY("G_py", x+1);
|
---|
149 | for (int y=0; y<ny; y++)
|
---|
150 | {
|
---|
151 | const Float_t s = gx->Integral(1, y+1);
|
---|
152 | const Float_t b = hx->Integral(1, y+1);
|
---|
153 | const Float_t sig = MMath::SignificanceLiMa(s+b, b);
|
---|
154 | res2.SetBinContent(x+1, y+1, sig);
|
---|
155 | }
|
---|
156 | delete hx;
|
---|
157 | delete gx;
|
---|
158 | }
|
---|
159 |
|
---|
160 | TGraph gr3;
|
---|
161 | TGraph gr4;
|
---|
162 | gr4.SetTitle("Significance Li/Ma vs. Hadronness-cut");
|
---|
163 |
|
---|
164 | TH1 *hx = h.ProjectionY("H_py");
|
---|
165 | TH1 *gx = g.ProjectionY("G_py");
|
---|
166 | for (int y=0; y<ny; y++)
|
---|
167 | {
|
---|
168 | const Float_t s = gx->Integral(1, y+1);
|
---|
169 | const Float_t b = hx->Integral(1, y+1);
|
---|
170 | const Float_t sig1 = MMath::SignificanceLiMa(s+b, b);
|
---|
171 | const Float_t sig2 = s<1 ? 0 : MMath::SignificanceLiMa(s+b, b)*TMath::Log10(s);
|
---|
172 |
|
---|
173 | gr3.SetPoint(y, h.GetYaxis()->GetBinLowEdge(y+2), sig1);
|
---|
174 | gr4.SetPoint(y, h.GetYaxis()->GetBinLowEdge(y+2), sig2);
|
---|
175 | }
|
---|
176 | delete hx;
|
---|
177 | delete gx;
|
---|
178 |
|
---|
179 | TCanvas &c = fDisplay->AddTab("OptCut");
|
---|
180 | c.SetBorderMode(0);
|
---|
181 | c.Divide(2,2);
|
---|
182 |
|
---|
183 | c.cd(1);
|
---|
184 | gPad->SetBorderMode(0);
|
---|
185 | gPad->SetFrameBorderMode(0);
|
---|
186 | gPad->SetLogx();
|
---|
187 | gPad->SetGridx();
|
---|
188 | gPad->SetGridy();
|
---|
189 | h.DrawCopy();
|
---|
190 | g.DrawCopy("same");
|
---|
191 | gr1.SetMarkerStyle(kFullDotMedium);
|
---|
192 | gr1.DrawClone("LP")->SetBit(kCanDelete);
|
---|
193 | gr2.SetLineColor(kBlue);
|
---|
194 | gr2.SetMarkerStyle(kFullDotMedium);
|
---|
195 | gr2.DrawClone("LP")->SetBit(kCanDelete);
|
---|
196 |
|
---|
197 | c.cd(3);
|
---|
198 | gPad->SetBorderMode(0);
|
---|
199 | gPad->SetFrameBorderMode(0);
|
---|
200 | gPad->SetGridx();
|
---|
201 | gPad->SetGridy();
|
---|
202 | gr4.SetMarkerStyle(kFullDotMedium);
|
---|
203 | gr4.DrawClone("ALP")->SetBit(kCanDelete);
|
---|
204 | gr3.SetLineColor(kBlue);
|
---|
205 | gr3.SetMarkerStyle(kFullDotMedium);
|
---|
206 | gr3.DrawClone("LP")->SetBit(kCanDelete);
|
---|
207 |
|
---|
208 | c.cd(2);
|
---|
209 | gPad->SetBorderMode(0);
|
---|
210 | gPad->SetFrameBorderMode(0);
|
---|
211 | gPad->SetLogx();
|
---|
212 | gPad->SetGridx();
|
---|
213 | gPad->SetGridy();
|
---|
214 | gPad->AddExec("color", "gStyle->SetPalette(1, 0);");
|
---|
215 | res1.SetMaximum(7);
|
---|
216 | res1.DrawCopy("colz");
|
---|
217 |
|
---|
218 | c.cd(4);
|
---|
219 | gPad->SetBorderMode(0);
|
---|
220 | gPad->SetFrameBorderMode(0);
|
---|
221 | gPad->SetLogx();
|
---|
222 | gPad->SetGridx();
|
---|
223 | gPad->SetGridy();
|
---|
224 | gPad->AddExec("color", "gStyle->SetPalette(1, 0);");
|
---|
225 | res2.SetMaximum(res2.GetMaximum()*1.075);
|
---|
226 | res2.DrawCopy("colz");
|
---|
227 | }
|
---|
228 |
|
---|
229 | /*
|
---|
230 | Bool_t MJSpectrum::InitWeighting(const MDataSet &set, MMcSpectrumWeight &w) const
|
---|
231 | {
|
---|
232 | fLog->Separator("Initialize energy weighting");
|
---|
233 |
|
---|
234 | if (!CheckEnv(w))
|
---|
235 | {
|
---|
236 | *fLog << err << "ERROR - Reading resources for MMcSpectrumWeight failed." << endl;
|
---|
237 | return kFALSE;
|
---|
238 | }
|
---|
239 |
|
---|
240 | TChain chain("RunHeaders");
|
---|
241 | set.AddFilesOn(chain);
|
---|
242 |
|
---|
243 | MMcCorsikaRunHeader *h=0;
|
---|
244 | chain.SetBranchAddress("MMcCorsikaRunHeader.", &h);
|
---|
245 | chain.GetEntry(1);
|
---|
246 |
|
---|
247 | if (!h)
|
---|
248 | {
|
---|
249 | *fLog << err << "ERROR - Couldn't read MMcCorsikaRunHeader from DataSet." << endl;
|
---|
250 | return kFALSE;
|
---|
251 | }
|
---|
252 |
|
---|
253 | if (!w.Set(*h))
|
---|
254 | {
|
---|
255 | *fLog << err << "ERROR - Initializing MMcSpectrumWeight failed." << endl;
|
---|
256 | return kFALSE;
|
---|
257 | }
|
---|
258 |
|
---|
259 | w.Print();
|
---|
260 | return kTRUE;
|
---|
261 | }
|
---|
262 |
|
---|
263 | Bool_t MJSpectrum::ReadOrigMCDistribution(const MDataSet &set, TH1 &h, MMcSpectrumWeight &weight) const
|
---|
264 | {
|
---|
265 | // Some debug output
|
---|
266 | fLog->Separator("Compiling original MC distribution");
|
---|
267 |
|
---|
268 | weight.SetNameMcEvt("MMcEvtBasic");
|
---|
269 | const TString w(weight.GetFormulaWeights());
|
---|
270 | weight.SetNameMcEvt();
|
---|
271 |
|
---|
272 | *fLog << inf << "Using weights: " << w << endl;
|
---|
273 | *fLog << "Please stand by, this may take a while..." << flush;
|
---|
274 |
|
---|
275 | if (fDisplay)
|
---|
276 | fDisplay->SetStatusLine1("Compiling MC distribution...");
|
---|
277 |
|
---|
278 | // Create chain
|
---|
279 | TChain chain("OriginalMC");
|
---|
280 | set.AddFilesOn(chain);
|
---|
281 |
|
---|
282 | // Prepare histogram
|
---|
283 | h.Reset();
|
---|
284 |
|
---|
285 | // Fill histogram from chain
|
---|
286 | h.SetDirectory(gROOT);
|
---|
287 | if (h.InheritsFrom(TH2::Class()))
|
---|
288 | {
|
---|
289 | h.SetNameTitle("ThetaEMC", "Event-Distribution vs Theta and Energy for MC (produced)");
|
---|
290 | h.SetXTitle("\\Theta [\\circ]");
|
---|
291 | h.SetYTitle("E [GeV]");
|
---|
292 | h.SetZTitle("Counts");
|
---|
293 | chain.Draw("MMcEvtBasic.fEnergy:MMcEvtBasic.fTelescopeTheta*TMath::RadToDeg()>>ThetaEMC", w, "goff");
|
---|
294 | }
|
---|
295 | else
|
---|
296 | {
|
---|
297 | h.SetNameTitle("ThetaMC", "Event-Distribution vs Theta for MC (produced)");
|
---|
298 | h.SetXTitle("\\Theta [\\circ]");
|
---|
299 | h.SetYTitle("Counts");
|
---|
300 | chain.Draw("MMcEvtBasic.fTelescopeTheta*TMath::RadToDeg()>>ThetaMC", w, "goff");
|
---|
301 | }
|
---|
302 | h.SetDirectory(0);
|
---|
303 |
|
---|
304 | *fLog << "done." << endl;
|
---|
305 | if (fDisplay)
|
---|
306 | fDisplay->SetStatusLine2("done.");
|
---|
307 |
|
---|
308 | if (h.GetEntries()>0)
|
---|
309 | return kTRUE;
|
---|
310 |
|
---|
311 | *fLog << err << "ERROR - Histogram with original MC distribution empty..." << endl;
|
---|
312 |
|
---|
313 | return h.GetEntries()>0;
|
---|
314 | }
|
---|
315 | */
|
---|
316 |
|
---|
317 | Bool_t MJTrainSeparation::GetEventsProduced(MDataSet &set, Double_t &num, Double_t &min, Double_t &max) const
|
---|
318 | {
|
---|
319 | TChain chain("OriginalMC");
|
---|
320 | set.AddFilesOn(chain);
|
---|
321 |
|
---|
322 | min = chain.GetMinimum("MMcEvtBasic.fEnergy");
|
---|
323 | max = chain.GetMaximum("MMcEvtBasic.fEnergy");
|
---|
324 |
|
---|
325 | num = chain.GetEntries();
|
---|
326 |
|
---|
327 | if (num<100)
|
---|
328 | *fLog << err << "ERROR - Less than 100 entries in OriginalMC-Tree of MC-Train-Data found." << endl;
|
---|
329 |
|
---|
330 | return num>=100;
|
---|
331 | }
|
---|
332 |
|
---|
333 | Double_t MJTrainSeparation::GetDataRate(MDataSet &set, Double_t &num) const
|
---|
334 | {
|
---|
335 | TChain chain1("Events");
|
---|
336 | set.AddFilesOff(chain1);
|
---|
337 |
|
---|
338 | num = chain1.GetEntries();
|
---|
339 | if (num<100)
|
---|
340 | {
|
---|
341 | *fLog << err << "ERROR - Less than 100 entries in Events-Tree of Train-Data found." << endl;
|
---|
342 | return -1;
|
---|
343 | }
|
---|
344 |
|
---|
345 | TChain chain("EffectiveOnTime");
|
---|
346 | set.AddFilesOff(chain);
|
---|
347 |
|
---|
348 | chain.Draw("MEffectiveOnTime.fVal", "MEffectiveOnTime.fVal", "goff");
|
---|
349 |
|
---|
350 | TH1 *h = dynamic_cast<TH1*>(gROOT->FindObject("htemp"));
|
---|
351 | if (!h)
|
---|
352 | {
|
---|
353 | *fLog << err << "ERROR - Weird things are happening (htemp not found)!" << endl;
|
---|
354 | return -1;
|
---|
355 | }
|
---|
356 |
|
---|
357 | const Double_t ontime = h->Integral();
|
---|
358 | delete h;
|
---|
359 |
|
---|
360 | if (ontime<1)
|
---|
361 | {
|
---|
362 | *fLog << err << "ERROR - Less than 1s of effective observation time found in Train-Data." << endl;
|
---|
363 | return -1;
|
---|
364 | }
|
---|
365 |
|
---|
366 | return num/ontime;
|
---|
367 | }
|
---|
368 |
|
---|
369 | Double_t MJTrainSeparation::GetNumMC(MDataSet &set) const
|
---|
370 | {
|
---|
371 | TChain chain1("Events");
|
---|
372 | set.AddFilesOn(chain1);
|
---|
373 |
|
---|
374 | const Double_t num = chain1.GetEntries();
|
---|
375 | if (num<100)
|
---|
376 | {
|
---|
377 | *fLog << err << "ERROR - Less than 100 entries in Events-Tree of Train-Data found." << endl;
|
---|
378 | return -1;
|
---|
379 | }
|
---|
380 |
|
---|
381 | return num;
|
---|
382 | }
|
---|
383 |
|
---|
384 | Bool_t MJTrainSeparation::AutoTrain(MDataSet &set, UInt_t &seton, UInt_t &setoff)
|
---|
385 | {
|
---|
386 | Double_t num, min, max;
|
---|
387 | if (!GetEventsProduced(set, num, min, max))
|
---|
388 | return kFALSE;
|
---|
389 |
|
---|
390 | *fLog << inf << "Using build-in radius of 300m to calculate collection area!" << endl;
|
---|
391 |
|
---|
392 | // Target spectrum
|
---|
393 | TF1 flx("Flux", "[0]/1000*(x/1000)^(-2.6)", min, max);
|
---|
394 | flx.SetParameter(0, 2e-7);
|
---|
395 |
|
---|
396 | // Number n0 of events this spectrum would produce per s and m^2
|
---|
397 | const Double_t n0 = flx.Integral(min, max); //[#]
|
---|
398 |
|
---|
399 | // Area produced in MC
|
---|
400 | const Double_t A = TMath::Pi()*300*300; //[m²]
|
---|
401 |
|
---|
402 | // Rate R of events this spectrum would produce per s
|
---|
403 | const Double_t R = n0*A; //[Hz]
|
---|
404 |
|
---|
405 | *fLog << "Gamma rate from the source inside the MC production area: " << R << "Hz" << endl;
|
---|
406 |
|
---|
407 | // Number N of events produced (in trainings sample)
|
---|
408 | const Double_t N = num; //[#]
|
---|
409 |
|
---|
410 | *fLog << "Events produced by MC inside the production area: " << TMath::Nint(num) << endl;
|
---|
411 |
|
---|
412 | // This correponds to an observation time T [s]
|
---|
413 | const Double_t T = N/R; //[s]
|
---|
414 |
|
---|
415 | *fLog << "Total time produced by the Monte Carlo: " << T << "s" << endl;
|
---|
416 |
|
---|
417 | // With an average data rate after star of
|
---|
418 | Double_t data=0;
|
---|
419 | const Double_t r = GetDataRate(set, data); //[Hz]
|
---|
420 |
|
---|
421 | *fLog << "Events measured per second effective on time: " << r << "Hz" << endl;
|
---|
422 | *fLog << "Total effective on time: " << data/r << "s" << endl;
|
---|
423 |
|
---|
424 | // this yields a number of n events to be read for training
|
---|
425 | const Double_t n = r*T; //[#]
|
---|
426 |
|
---|
427 | *fLog << "Events to be read from the data sample: " << TMath::Nint(n) << endl;
|
---|
428 | *fLog << "Events available in data sample: " << data << endl;
|
---|
429 |
|
---|
430 | if (r<0)
|
---|
431 | return kFALSE;
|
---|
432 |
|
---|
433 | Double_t nummc = GetNumMC(set);
|
---|
434 |
|
---|
435 | *fLog << "Events available in MC sample: " << nummc << endl;
|
---|
436 |
|
---|
437 | *fLog << "MC read probability: " << data/n << endl;
|
---|
438 |
|
---|
439 | // more data requested than available => Scale down num MC events
|
---|
440 | Double_t on, off;
|
---|
441 | if (data<n)
|
---|
442 | {
|
---|
443 | on = TMath::Nint(nummc*data/n);
|
---|
444 | off = TMath::Nint(data);
|
---|
445 | *fLog << warn;
|
---|
446 | *fLog << "Not enough data events available... scaling by " << data/n << endl;
|
---|
447 | *fLog << inf;
|
---|
448 | }
|
---|
449 | else
|
---|
450 | {
|
---|
451 | on = TMath::Nint(nummc);
|
---|
452 | off = TMath::Nint(n);
|
---|
453 | }
|
---|
454 |
|
---|
455 | if (seton>0 && seton<on)
|
---|
456 | {
|
---|
457 | setoff = TMath::Nint(off*seton/on);
|
---|
458 | *fLog << "Less MC events requested... scaling by " << seton/on << endl;
|
---|
459 | }
|
---|
460 | else
|
---|
461 | {
|
---|
462 | seton = TMath::Nint(on);
|
---|
463 | setoff = TMath::Nint(off);
|
---|
464 | }
|
---|
465 |
|
---|
466 | *fLog << "Target number of MC events: " << seton << endl;
|
---|
467 | *fLog << "Target number of data events: " << setoff << endl;
|
---|
468 |
|
---|
469 | /*
|
---|
470 | An event rate dependent selection?
|
---|
471 | ----------------------------------
|
---|
472 | Total average data rate: R
|
---|
473 | Goal number of events: N
|
---|
474 | Number of data events: N0
|
---|
475 | Rate assigned to single evt: r
|
---|
476 |
|
---|
477 | Selection probability: N/N0 * r/R
|
---|
478 |
|
---|
479 | f := N/N0 * r
|
---|
480 |
|
---|
481 | MF f("f * MEventRate.fRate < rand");
|
---|
482 | */
|
---|
483 |
|
---|
484 | return kTRUE;
|
---|
485 | }
|
---|
486 |
|
---|
487 | Bool_t MJTrainSeparation::Train(const char *out)
|
---|
488 | {
|
---|
489 | if (!fDataSetTrain.IsValid())
|
---|
490 | {
|
---|
491 | *fLog << err << "ERROR - DataSet for training invalid!" << endl;
|
---|
492 | return kFALSE;
|
---|
493 | }
|
---|
494 | if (!fDataSetTest.IsValid())
|
---|
495 | {
|
---|
496 | *fLog << err << "ERROR - DataSet for testing invalid!" << endl;
|
---|
497 | return kFALSE;
|
---|
498 | }
|
---|
499 |
|
---|
500 | // ----------------------- Auto Train? ----------------------
|
---|
501 |
|
---|
502 | if (fAutoTrain)
|
---|
503 | {
|
---|
504 | fLog->Separator("Auto-Training -- Train-Data");
|
---|
505 | if (!AutoTrain(fDataSetTrain, fNumTrainOn, fNumTrainOff))
|
---|
506 | return kFALSE;
|
---|
507 | fLog->Separator("Auto-Training -- Test-Data");
|
---|
508 | if (!AutoTrain(fDataSetTest, fNumTestOn, fNumTestOff))
|
---|
509 | return kFALSE;
|
---|
510 | }
|
---|
511 |
|
---|
512 | // --------------------- Setup files --------------------
|
---|
513 | MReadMarsFile read1("Events");
|
---|
514 | MReadMarsFile read2("Events");
|
---|
515 | MReadMarsFile read3("Events");
|
---|
516 | MReadMarsFile read4("Events");
|
---|
517 | read1.DisableAutoScheme();
|
---|
518 | read2.DisableAutoScheme();
|
---|
519 | read3.DisableAutoScheme();
|
---|
520 | read4.DisableAutoScheme();
|
---|
521 |
|
---|
522 | fDataSetTrain.AddFilesOn(read1);
|
---|
523 | fDataSetTrain.AddFilesOff(read3);
|
---|
524 |
|
---|
525 | fDataSetTest.AddFilesOff(read2);
|
---|
526 | fDataSetTest.AddFilesOn(read4);
|
---|
527 |
|
---|
528 | // ----------------------- Setup RF ----------------------
|
---|
529 | MHMatrix train("Train");
|
---|
530 | train.AddColumns(fRules);
|
---|
531 | train.AddColumn("MHadronness.fVal");
|
---|
532 |
|
---|
533 | // ----------------------- Fill Matrix RF ----------------------
|
---|
534 |
|
---|
535 | MParameterD had("MHadronness");
|
---|
536 |
|
---|
537 | MParList plistx;
|
---|
538 | plistx.AddToList(&had);
|
---|
539 | plistx.AddToList(this);
|
---|
540 |
|
---|
541 | MTFillMatrix fill;
|
---|
542 | fill.SetLogStream(fLog);
|
---|
543 | fill.SetDisplay(fDisplay);
|
---|
544 | fill.AddPreCuts(fPreCuts);
|
---|
545 | fill.AddPreCuts(fTrainCuts);
|
---|
546 |
|
---|
547 | // Set classifier for gammas
|
---|
548 | had.SetVal(0);
|
---|
549 | fill.SetName("FillGammas");
|
---|
550 | fill.SetDestMatrix1(&train, fNumTrainOn);
|
---|
551 | fill.SetReader(&read1);
|
---|
552 | if (!fill.Process(plistx))
|
---|
553 | return kFALSE;
|
---|
554 |
|
---|
555 | const Int_t numgammastrn = train.GetNumRows();
|
---|
556 | if (numgammastrn==0)
|
---|
557 | {
|
---|
558 | *fLog << err << "ERROR - No gammas available for training... aborting." << endl;
|
---|
559 | return kFALSE;
|
---|
560 | }
|
---|
561 |
|
---|
562 | // Set classifier for hadrons
|
---|
563 | had.SetVal(1);
|
---|
564 | fill.SetName("FillBackground");
|
---|
565 | fill.SetDestMatrix1(&train, fNumTrainOff);
|
---|
566 | fill.SetReader(&read3);
|
---|
567 | if (!fill.Process(plistx))
|
---|
568 | return kFALSE;
|
---|
569 |
|
---|
570 | const Int_t numbackgrndtrn = train.GetNumRows()-numgammastrn;
|
---|
571 | if (numbackgrndtrn==0)
|
---|
572 | {
|
---|
573 | *fLog << err << "ERROR - No background available for training... aborting." << endl;
|
---|
574 | return kFALSE;
|
---|
575 | }
|
---|
576 |
|
---|
577 | // ------------------------ Train RF --------------------------
|
---|
578 |
|
---|
579 | MRanForestCalc rf;
|
---|
580 | rf.SetNumTrees(fNumTrees);
|
---|
581 | rf.SetNdSize(fNdSize);
|
---|
582 | rf.SetNumTry(fNumTry);
|
---|
583 | rf.SetNumObsoleteVariables(1);
|
---|
584 | rf.SetDebug(fDebug);
|
---|
585 | rf.SetDisplay(fDisplay);
|
---|
586 | rf.SetLogStream(fLog);
|
---|
587 | rf.SetFileName(out);
|
---|
588 | rf.SetNameOutput("MHadronness");
|
---|
589 |
|
---|
590 | if (fUseRegression)
|
---|
591 | {
|
---|
592 | if (!rf.TrainSingleRF(train)) // regression
|
---|
593 | return kFALSE;
|
---|
594 | }
|
---|
595 | else
|
---|
596 | {
|
---|
597 | MBinning b(2, -0.5, 1.5, "BinningHadronness", "lin");
|
---|
598 | if (!rf.TrainSingleRF(train, b.GetEdgesD())) // classification
|
---|
599 | return kFALSE;
|
---|
600 | }
|
---|
601 |
|
---|
602 | //if (!rf.TrainMultiRF(train, b.GetEdgesD())) // classification
|
---|
603 | // return;
|
---|
604 |
|
---|
605 | //fDisplay = rf.GetDisplay();
|
---|
606 |
|
---|
607 |
|
---|
608 | *fLog << all;
|
---|
609 | fLog->Separator("The forest was tested with...");
|
---|
610 |
|
---|
611 | *fLog << "Training method:" << endl;
|
---|
612 | *fLog << " * " << (fUseRegression?"regression":"classification") << endl;
|
---|
613 | *fLog << endl;
|
---|
614 | *fLog << "Events used for training:" << endl;
|
---|
615 | *fLog << " * Gammas: " << numgammastrn << endl;
|
---|
616 | *fLog << " * Background: " << numbackgrndtrn << endl;
|
---|
617 | *fLog << endl;
|
---|
618 | *fLog << "Gamma/Background ratio:" << endl;
|
---|
619 | *fLog << " * Requested: " << (float)fNumTrainOn/fNumTrainOff << endl;
|
---|
620 | *fLog << " * Result: " << (float)numgammastrn/numbackgrndtrn << endl;
|
---|
621 |
|
---|
622 | if (!fDataSetTest.IsValid())
|
---|
623 | return kTRUE;
|
---|
624 |
|
---|
625 | // --------------------- Display result ----------------------
|
---|
626 | fLog->Separator("Test");
|
---|
627 |
|
---|
628 | MParList plist;
|
---|
629 | MTaskList tlist;
|
---|
630 | plist.AddToList(this);
|
---|
631 | plist.AddToList(&tlist);
|
---|
632 |
|
---|
633 | MMcEvt mcevt;
|
---|
634 | plist.AddToList(&mcevt);
|
---|
635 |
|
---|
636 | // ----- Setup histograms -----
|
---|
637 | MBinning binsy(50, 0 , 1, "BinningMH3Y", "lin");
|
---|
638 | MBinning binsx(40, 10, 100000, "BinningMH3X", "log");
|
---|
639 |
|
---|
640 | plist.AddToList(&binsx);
|
---|
641 | plist.AddToList(&binsy);
|
---|
642 |
|
---|
643 | MH3 h31("MHillas.fSize", "MHadronness.fVal");
|
---|
644 | MH3 h32("MHillas.fSize", "MHadronness.fVal");
|
---|
645 | MH3 h40("MMcEvt.fEnergy", "MHadronness.fVal");
|
---|
646 | h31.SetTitle("Background probability vs. Size:Size [phe]:Hadronness");
|
---|
647 | h32.SetTitle("Background probability vs. Size:Size [phe]:Hadronness");
|
---|
648 | h40.SetTitle("Background probability vs. Energy:Energy [GeV]:Hadronness");
|
---|
649 |
|
---|
650 | MHHadronness hist;
|
---|
651 |
|
---|
652 | // ----- Setup tasks -----
|
---|
653 | MFillH fillh0(&hist, "", "FillHadronness");
|
---|
654 | MFillH fillh1(&h31);
|
---|
655 | MFillH fillh2(&h32);
|
---|
656 | MFillH fillh4(&h40);
|
---|
657 | fillh1.SetNameTab("Background");
|
---|
658 | fillh2.SetNameTab("GammasH");
|
---|
659 | fillh4.SetNameTab("GammasE");
|
---|
660 | fillh0.SetBit(MFillH::kDoNotDisplay);
|
---|
661 |
|
---|
662 | // ----- Setup filter -----
|
---|
663 | MFilterList precuts;
|
---|
664 | precuts.AddToList(fPreCuts);
|
---|
665 | precuts.AddToList(fTestCuts);
|
---|
666 |
|
---|
667 | MContinue c0(&precuts);
|
---|
668 | c0.SetName("PreCuts");
|
---|
669 | c0.SetInverted();
|
---|
670 |
|
---|
671 | MFEventSelector sel;
|
---|
672 | sel.SetNumSelectEvts(fNumTestOff);
|
---|
673 |
|
---|
674 | MContinue c1(&sel);
|
---|
675 | c1.SetInverted();
|
---|
676 |
|
---|
677 | // ----- Setup tasklist -----
|
---|
678 | tlist.AddToList(&read2);
|
---|
679 | tlist.AddToList(&c0);
|
---|
680 | tlist.AddToList(&c1);
|
---|
681 | tlist.AddToList(&rf);
|
---|
682 | tlist.AddToList(&fillh0);
|
---|
683 | tlist.AddToList(&fillh1);
|
---|
684 |
|
---|
685 | // ----- Run eventloop on gammas -----
|
---|
686 | MEvtLoop loop;
|
---|
687 | loop.SetDisplay(fDisplay);
|
---|
688 | loop.SetLogStream(fLog);
|
---|
689 | loop.SetParList(&plist);
|
---|
690 |
|
---|
691 | if (!loop.Eventloop())
|
---|
692 | return kFALSE;
|
---|
693 |
|
---|
694 | // ----- Setup and run eventloop on background -----
|
---|
695 | sel.SetNumSelectEvts(fNumTestOn);
|
---|
696 | fillh0.ResetBit(MFillH::kDoNotDisplay);
|
---|
697 |
|
---|
698 | tlist.Replace(&read4);
|
---|
699 | tlist.Replace(&fillh2);
|
---|
700 | tlist.AddToListAfter(&fillh4, &fillh2);
|
---|
701 |
|
---|
702 | if (!loop.Eventloop())
|
---|
703 | return kFALSE;
|
---|
704 |
|
---|
705 | DisplayResult(h31, h32);
|
---|
706 |
|
---|
707 | if (!WriteDisplay(out))
|
---|
708 | return kFALSE;
|
---|
709 |
|
---|
710 | *fLog << all;
|
---|
711 | fLog->Separator("The forest was tested with...");
|
---|
712 |
|
---|
713 | const Double_t numgammastst = h32.GetHist().GetEntries();
|
---|
714 | const Double_t numbackgrndtst = h31.GetHist().GetEntries();
|
---|
715 |
|
---|
716 | *fLog << "Events used for test:" << endl;
|
---|
717 | *fLog << " * Gammas: " << numgammastst << endl;
|
---|
718 | *fLog << " * Background: " << numbackgrndtst << endl;
|
---|
719 | *fLog << endl;
|
---|
720 | *fLog << "Gamma/Background ratio:" << endl;
|
---|
721 | *fLog << " * Requested: " << (float)fNumTestOn/fNumTestOff << endl;
|
---|
722 | *fLog << " * Result: " << (float)numgammastst/numbackgrndtst << endl;
|
---|
723 |
|
---|
724 | return kTRUE;
|
---|
725 | }
|
---|
726 |
|
---|