Context Navigation

source: trunk/MagicSoft/Mars/mjtrain/MJTrainSeparation.cc@ 8154

Visit:

Last change on this file since 8154 was 7760, checked in by tbretz, 18 years ago
* empty log message *
File size: 27.0 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, 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 <TMarker.h>
37	#include <TCanvas.h>
38	#include <TStopwatch.h>
39	#include <TVirtualPad.h>
40
41	#include "MHMatrix.h"
42
43	#include "MLog.h"
44	#include "MLogManip.h"
45
46	// tools
47	#include "MMath.h"
48	#include "MDataSet.h"
49	#include "MTFillMatrix.h"
50	#include "MStatusDisplay.h"
51
52	// eventloop
53	#include "MParList.h"
54	#include "MTaskList.h"
55	#include "MEvtLoop.h"
56
57	// tasks
58	#include "MReadMarsFile.h"
59	#include "MContinue.h"
60	#include "MFillH.h"
61	#include "MSrcPosRndm.h"
62	#include "MHillasCalc.h"
63	#include "MRanForestCalc.h"
64	#include "MParameterCalc.h"
65
66	// container
67	#include "MMcEvt.hxx"
68	#include "MParameters.h"
69
70	// histograms
71	#include "MBinning.h"
72	#include "MH3.h"
73	#include "MHHadronness.h"
74
75	// filter
76	#include "MF.h"
77	#include "MFEventSelector.h"
78	#include "MFilterList.h"
79
80	ClassImp(MJTrainSeparation);
81
82	using namespace std;
83
84	void MJTrainSeparation::DisplayResult(MH3 &h31, MH3 &h32, Float_t ontime)
85	{
86	TH2D &g = (TH2D&)h32.GetHist();
87	TH2D &h = (TH2D&)h31.GetHist();
88
89	h.SetMarkerColor(kRed);
90	g.SetMarkerColor(kGreen);
91
92	TH2D res1(g);
93	TH2D res2(g);
94
95	h.SetTitle("Hadronness-Distribution vs. Size");
96	res1.SetTitle("Significance Li/Ma");
97	res1.SetXTitle("Size [phe]");
98	res1.SetYTitle("Hadronness");
99	res2.SetTitle("Significance-Distribution");
100	res2.SetXTitle("Size-Cut [phe]");
101	res2.SetYTitle("Hadronness-Cut");
102	res1.SetContour(50);
103	res2.SetContour(50);
104
105	const Int_t nx = h.GetNbinsX();
106	const Int_t ny = h.GetNbinsY();
107
108	gROOT->SetSelectedPad(NULL);
109
110
111	Double_t Stot = 0;
112	Double_t Btot = 0;
113
114	Double_t max2 = -1;
115
116	TGraph gr1;
117	TGraph gr2;
118	for (int x=nx-1; x>=0; x--)
119	{
120	TH1 *hx = h.ProjectionY("H_py", x+1, x+1);
121	TH1 *gx = g.ProjectionY("G_py", x+1, x+1);
122
123	Double_t S = 0;
124	Double_t B = 0;
125
126	Double_t max1 = -1;
127	Int_t maxy1 = 0;
128	Int_t maxy2 = 0;
129	for (int y=ny-1; y>=0; y--)
130	{
131	const Float_t s = gx->Integral(1, y+1);
132	const Float_t b = hx->Integral(1, y+1);
133	const Float_t sig1 = MMath::SignificanceLiMa(s+b, b);
134	const Float_t sig2 = MMath::SignificanceLiMa(s+Stot+b+Btot, b+Btot)*TMath::Log10(s+Stot+1);
135	if (sig1>max1)
136	{
137	maxy1 = y;
138	max1 = sig1;
139	}
140	if (sig2>max2)
141	{
142	maxy2 = y;
143	max2 = sig2;
144
145	S=s;
146	B=b;
147	}
148
149	res1.SetBinContent(x+1, y+1, sig1);
150	}
151
152	Stot += S;
153	Btot += B;
154
155	gr1.SetPoint(x, h.GetXaxis()->GetBinCenter(x+1), h.GetYaxis()->GetBinCenter(maxy1+1));
156	gr2.SetPoint(x, h.GetXaxis()->GetBinCenter(x+1), h.GetYaxis()->GetBinCenter(maxy2+1));
157
158	delete hx;
159	delete gx;
160	}
161
162	//cout << "--> " << MMath::SignificanceLiMa(Stot+Btot, Btot) << " ";
163	//cout << Stot << " " << Btot << endl;
164
165
166	Int_t mx1=0;
167	Int_t my1=0;
168	Int_t mx2=0;
169	Int_t my2=0;
170	Int_t s1=0;
171	Int_t b1=0;
172	Int_t s2=0;
173	Int_t b2=0;
174	Double_t sig1=-1;
175	Double_t sig2=-1;
176	for (int x=0; x<nx; x++)
177	{
178	TH1 *hx = h.ProjectionY("H_py", x+1);
179	TH1 *gx = g.ProjectionY("G_py", x+1);
180	for (int y=0; y<ny; y++)
181	{
182	const Float_t s = gx->Integral(1, y+1);
183	const Float_t b = hx->Integral(1, y+1);
184	const Float_t sig = MMath::SignificanceLiMa(s+b, b);
185	res2.SetBinContent(x+1, y+1, sig);
186
187	// Search for top-rightmost maximum
188	if (sig>=sig1)
189	{
190	mx1=x+1;
191	my1=y+1;
192	s1 = TMath::Nint(s);
193	b1 = TMath::Nint(b);
194	sig1=sig;
195	}
196	if (TMath::Log10(s)*sig>=sig2)
197	{
198	mx2=x+1;
199	my2=y+1;
200	s2 = TMath::Nint(s);
201	b2 = TMath::Nint(b);
202	sig2=TMath::Log10(s)*sig;
203	}
204	}
205	delete hx;
206	delete gx;
207	}
208
209	TGraph gr3;
210	TGraph gr4;
211	gr4.SetTitle("Significance Li/Ma vs. Hadronness-cut");
212
213	TH1 *hx = h.ProjectionY("H_py");
214	TH1 *gx = g.ProjectionY("G_py");
215	for (int y=0; y<ny; y++)
216	{
217	const Float_t s = gx->Integral(1, y+1);
218	const Float_t b = hx->Integral(1, y+1);
219	const Float_t sig1 = MMath::SignificanceLiMa(s+b, b);
220	const Float_t sig2 = s<1 ? 0 : MMath::SignificanceLiMa(s+b, b)*TMath::Log10(s);
221
222	gr3.SetPoint(y, h.GetYaxis()->GetBinLowEdge(y+2), sig1);
223	gr4.SetPoint(y, h.GetYaxis()->GetBinLowEdge(y+2), sig2);
224	}
225	delete hx;
226	delete gx;
227
228	if (fDisplay)
229	{
230	TCanvas &c = fDisplay->AddTab("OptCut");
231	c.SetBorderMode(0);
232	c.Divide(2,2);
233
234	c.cd(1);
235	gPad->SetBorderMode(0);
236	gPad->SetFrameBorderMode(0);
237	gPad->SetLogx();
238	gPad->SetGridx();
239	gPad->SetGridy();
240	h.DrawCopy();
241	g.DrawCopy("same");
242	gr1.SetMarkerStyle(kFullDotMedium);
243	gr1.DrawClone("LP")->SetBit(kCanDelete);
244	gr2.SetLineColor(kBlue);
245	gr2.SetMarkerStyle(kFullDotMedium);
246	gr2.DrawClone("LP")->SetBit(kCanDelete);
247
248	c.cd(3);
249	gPad->SetBorderMode(0);
250	gPad->SetFrameBorderMode(0);
251	gPad->SetGridx();
252	gPad->SetGridy();
253	gr4.SetMinimum(0);
254	gr4.SetMarkerStyle(kFullDotMedium);
255	gr4.DrawClone("ALP")->SetBit(kCanDelete);
256	gr3.SetLineColor(kBlue);
257	gr3.SetMarkerStyle(kFullDotMedium);
258	gr3.DrawClone("LP")->SetBit(kCanDelete);
259
260	c.cd(2);
261	gPad->SetBorderMode(0);
262	gPad->SetFrameBorderMode(0);
263	gPad->SetLogx();
264	gPad->SetGridx();
265	gPad->SetGridy();
266	gPad->AddExec("color", "gStyle->SetPalette(1, 0);");
267	res1.SetMaximum(7);
268	res1.DrawCopy("colz");
269
270	c.cd(4);
271	gPad->SetBorderMode(0);
272	gPad->SetFrameBorderMode(0);
273	gPad->SetLogx();
274	gPad->SetGridx();
275	gPad->SetGridy();
276	gPad->AddExec("color", "gStyle->SetPalette(1, 0);");
277	res2.SetMaximum(res2.GetMaximum()*1.05);
278	res2.DrawCopy("colz");
279
280	// Int_t mx, my, mz;
281	// res2.GetMaximumBin(mx, my, mz);
282
283	TMarker m;
284	m.SetMarkerStyle(kStar);
285	m.DrawMarker(res2.GetXaxis()->GetBinCenter(mx1), res2.GetYaxis()->GetBinCenter(my1));
286	m.SetMarkerStyle(kPlus);
287	m.DrawMarker(res2.GetXaxis()->GetBinCenter(mx2), res2.GetYaxis()->GetBinCenter(my2));
288	}
289
290	*fLog << all << "Observation Time: " << TMath::Nint(ontime/60) << "min" << endl;
291	*fLog << "Maximum Significance: " << Form("%.1f", sig1) << " [";
292	*fLog << Form("%.1f", sig1/TMath::Sqrt(ontime/3600)) << "/sqrt(h)]";
293	*fLog << endl;
294
295	*fLog << "Significance: S=" << Form("%.1f", sig1) << " E=" << s1 << " B=" << b1 << " h<";
296	*fLog << Form("%.2f", res2.GetYaxis()->GetBinCenter(my1)) << " s>";
297	*fLog << Form("%3d", TMath::Nint(res2.GetXaxis()->GetBinCenter(mx1))) << endl;
298	fLog << "SignificanceLogE: S=" << Form("%.1f", sig2/TMath::Log10(s2)) << " E=" << s2 << " B=" << b2 << " h<";
299	*fLog << Form("%.2f", res2.GetYaxis()->GetBinCenter(my2)) << " s>";
300	*fLog << Form("%3d", TMath::Nint(res2.GetXaxis()->GetBinCenter(mx2))) << endl;
301	*fLog << endl;
302	}
303
304	/*
305	Bool_t MJSpectrum::InitWeighting(const MDataSet &set, MMcSpectrumWeight &w) const
306	{
307	fLog->Separator("Initialize energy weighting");
308
309	if (!CheckEnv(w))
310	{
311	*fLog << err << "ERROR - Reading resources for MMcSpectrumWeight failed." << endl;
312	return kFALSE;
313	}
314
315	TChain chain("RunHeaders");
316	set.AddFilesOn(chain);
317
318	MMcCorsikaRunHeader *h=0;
319	chain.SetBranchAddress("MMcCorsikaRunHeader.", &h);
320	chain.GetEntry(1);
321
322	if (!h)
323	{
324	*fLog << err << "ERROR - Couldn't read MMcCorsikaRunHeader from DataSet." << endl;
325	return kFALSE;
326	}
327
328	if (!w.Set(*h))
329	{
330	*fLog << err << "ERROR - Initializing MMcSpectrumWeight failed." << endl;
331	return kFALSE;
332	}
333
334	w.Print();
335	return kTRUE;
336	}
337
338	Bool_t MJSpectrum::ReadOrigMCDistribution(const MDataSet &set, TH1 &h, MMcSpectrumWeight &weight) const
339	{
340	// Some debug output
341	fLog->Separator("Compiling original MC distribution");
342
343	weight.SetNameMcEvt("MMcEvtBasic");
344	const TString w(weight.GetFormulaWeights());
345	weight.SetNameMcEvt();
346
347	*fLog << inf << "Using weights: " << w << endl;
348	*fLog << "Please stand by, this may take a while..." << flush;
349
350	if (fDisplay)
351	fDisplay->SetStatusLine1("Compiling MC distribution...");
352
353	// Create chain
354	TChain chain("OriginalMC");
355	set.AddFilesOn(chain);
356
357	// Prepare histogram
358	h.Reset();
359
360	// Fill histogram from chain
361	h.SetDirectory(gROOT);
362	if (h.InheritsFrom(TH2::Class()))
363	{
364	h.SetNameTitle("ThetaEMC", "Event-Distribution vs Theta and Energy for MC (produced)");
365	h.SetXTitle("\\Theta [\\circ]");
366	h.SetYTitle("E [GeV]");
367	h.SetZTitle("Counts");
368	chain.Draw("MMcEvtBasic.fEnergy:MMcEvtBasic.fTelescopeTheta*TMath::RadToDeg()>>ThetaEMC", w, "goff");
369	}
370	else
371	{
372	h.SetNameTitle("ThetaMC", "Event-Distribution vs Theta for MC (produced)");
373	h.SetXTitle("\\Theta [\\circ]");
374	h.SetYTitle("Counts");
375	chain.Draw("MMcEvtBasic.fTelescopeTheta*TMath::RadToDeg()>>ThetaMC", w, "goff");
376	}
377	h.SetDirectory(0);
378
379	*fLog << "done." << endl;
380	if (fDisplay)
381	fDisplay->SetStatusLine2("done.");
382
383	if (h.GetEntries()>0)
384	return kTRUE;
385
386	*fLog << err << "ERROR - Histogram with original MC distribution empty..." << endl;
387
388	return h.GetEntries()>0;
389	}
390	*/
391
392	Bool_t MJTrainSeparation::GetEventsProduced(MDataSet &set, Double_t &num, Double_t &min, Double_t &max) const
393	{
394	TChain chain("OriginalMC");
395	set.AddFilesOn(chain);
396
397	min = chain.GetMinimum("MMcEvtBasic.fEnergy");
398	max = chain.GetMaximum("MMcEvtBasic.fEnergy");
399
400	num = chain.GetEntries();
401
402	if (num<100)
403	*fLog << err << "ERROR - Less than 100 entries in OriginalMC-Tree of MC-Train-Data found." << endl;
404
405	return num>=100;
406	}
407
408	Double_t MJTrainSeparation::GetDataRate(MDataSet &set, Double_t &num) const
409	{
410	TChain chain1("Events");
411	set.AddFilesOff(chain1);
412
413	num = chain1.GetEntries();
414	if (num<100)
415	{
416	*fLog << err << "ERROR - Less than 100 entries in Events-Tree of Train-Data found." << endl;
417	return -1;
418	}
419
420	TChain chain("EffectiveOnTime");
421	set.AddFilesOff(chain);
422
423	chain.Draw("MEffectiveOnTime.fVal", "MEffectiveOnTime.fVal", "goff");
424
425	TH1 h = dynamic_cast<TH1>(gROOT->FindObject("htemp"));
426	if (!h)
427	{
428	*fLog << err << "ERROR - Weird things are happening (htemp not found)!" << endl;
429	return -1;
430	}
431
432	const Double_t ontime = h->Integral();
433	delete h;
434
435	if (ontime<1)
436	{
437	*fLog << err << "ERROR - Less than 1s of effective observation time found in Train-Data." << endl;
438	return -1;
439	}
440
441	return num/ontime;
442	}
443
444	Double_t MJTrainSeparation::GetNumMC(MDataSet &set) const
445	{
446	TChain chain1("Events");
447	set.AddFilesOn(chain1);
448
449	const Double_t num = chain1.GetEntries();
450	if (num<100)
451	{
452	*fLog << err << "ERROR - Less than 100 entries in Events-Tree of Train-Data found." << endl;
453	return -1;
454	}
455
456	return num;
457	}
458
459	Float_t MJTrainSeparation::AutoTrain(MDataSet &set, Type_t typon, Type_t typof, Float_t flux)
460	{
461	Double_t num, min, max;
462	if (!GetEventsProduced(set, num, min, max))
463	return -1;
464
465	*fLog << inf << "Using build-in radius of 300m to calculate collection area!" << endl;
466
467	// Target spectrum
468	TF1 flx("Flux", "[0]/1000*(x/1000)^(-2.6)", min, max);
469	flx.SetParameter(0, flux);
470
471	// Number n0 of events this spectrum would produce per s and m^2
472	const Double_t n0 = flx.Integral(min, max); //[#]
473
474	// Area produced in MC
475	const Double_t A = TMath::Pi()300300; //[m²]
476
477	// Rate R of events this spectrum would produce per s
478	const Double_t R = n0*A; //[Hz]
479
480	fLog << "Source Spectrum: " << flux << " (E/TeV)^(-2.6) * TeVm^2s" << endl;
481
482	*fLog << "Gamma rate from the source inside the MC production area: " << R << "Hz" << endl;
483
484	// Number N of events produced (in trainings sample)
485	const Double_t N = num; //[#]
486
487	*fLog << "Events produced by MC inside the production area: " << TMath::Nint(num) << endl;
488
489	// This correponds to an observation time T [s]
490	const Double_t T = N/R; //[s]
491
492	*fLog << "Total time produced by the Monte Carlo: " << T << "s" << endl;
493
494	// With an average data rate after star of
495	Double_t data=0;
496	const Double_t r = GetDataRate(set, data); //[Hz]
497	Double_t ontime = data/r;
498
499	*fLog << "Events measured per second effective on time: " << r << "Hz" << endl;
500	*fLog << "Total effective on time: " << ontime << "s" << endl;
501
502	const Double_t ratio = T/ontime;
503	*fLog << "Ratio of Monte Carlo to data observation time: " << ratio << endl;
504
505	// 3570.5/43440.2 = 0.082
506
507
508	// this yields a number of n events to be read for training
509	const Double_t n = r*T; //[#]
510
511	*fLog << "Events to be read from the data sample: " << TMath::Nint(n) << endl;
512	*fLog << "Events available in data sample: " << data << endl;
513
514	if (r<0)
515	return -1;
516
517	Double_t nummc = GetNumMC(set);
518
519	*fLog << "Events available in MC sample: " << nummc << endl;
520
521	// *fLog << "MC read probability: " << data/n << endl;
522
523	// more data requested than available => Scale down num MC events
524	Double_t on, off;
525	if (data<n)
526	{
527	on = TMath::Nint(nummc*data/n);
528	off = TMath::Nint(data);
529	*fLog << warn;
530	*fLog << "Not enough data events available... scaling MC to data by " << data/n << endl;
531	*fLog << inf;
532	}
533	else
534	{
535	on = TMath::Nint(nummc);
536	off = TMath::Nint(n);
537	}
538
539	if (fNum[typon]>0 && fNum[typon]<on)
540	{
541	fNum[typof] = TMath::Nint(off*fNum[typon]/on);
542	ontime *= fNum[typon]/on;
543	*fLog << warn << "Less MC events requested... scaling data to MC by " << fNum[typon]/on << endl;
544	}
545	else
546	{
547	fNum[typon] = TMath::Nint(on);
548	fNum[typof] = TMath::Nint(off);
549	}
550
551	*fLog << inf;
552	*fLog << "Target number of MC events: " << fNum[typon] << endl;
553	*fLog << "Target number of data events: " << fNum[typof] << endl;
554
555	/*
556	An event rate dependent selection?
557	----------------------------------
558	Total average data rate: R
559	Goal number of events: N
560	Number of data events: N0
561	Rate assigned to single evt: r
562
563	Selection probability: N/N0 * r/R
564
565	f := N/N0 * r
566
567	MF f("f * MEventRate.fRate < rand");
568	*/
569
570	return ontime;
571	}
572
573	Bool_t MJTrainSeparation::Train(const char *out)
574	{
575	if (!fDataSetTrain.IsValid())
576	{
577	*fLog << err << "ERROR - DataSet for training invalid!" << endl;
578	return kFALSE;
579	}
580	if (!fDataSetTest.IsValid())
581	{
582	*fLog << err << "ERROR - DataSet for testing invalid!" << endl;
583	return kFALSE;
584	}
585
586	if (fDataSetTrain.IsWobbleMode()!=fDataSetTest.IsWobbleMode())
587	{
588	*fLog << err << "ERROR - Train- and Test-DataSet have different observation modes!" << endl;
589	return kFALSE;
590	}
591
592	TStopwatch clock;
593	clock.Start();
594
595	// ----------------------- Auto Train? ----------------------
596
597	Float_t ontime = -1;
598	if (fAutoTrain)
599	{
600	fLog->Separator("Auto-Training -- Train-Data");
601	if (AutoTrain(fDataSetTrain, kTrainOn, kTrainOff, fFluxTrain)<0)
602	return kFALSE;
603	fLog->Separator("Auto-Training -- Test-Data");
604	ontime = AutoTrain(fDataSetTest, kTestOn, kTestOff, fFluxTest);
605	if (ontime<0)
606	return kFALSE;
607	}
608
609	// --------------------- Setup files --------------------
610	MReadMarsFile read1("Events");
611	MReadMarsFile read2("Events");
612	MReadMarsFile read3("Events");
613	MReadMarsFile read4("Events");
614	read1.DisableAutoScheme();
615	read2.DisableAutoScheme();
616	read3.DisableAutoScheme();
617	read4.DisableAutoScheme();
618
619	// Setup four reading tasks with the on- and off-data of the two datasets
620	fDataSetTrain.AddFilesOn(read1);
621	fDataSetTrain.AddFilesOff(read3);
622
623	fDataSetTest.AddFilesOff(read2);
624	fDataSetTest.AddFilesOn(read4);
625
626	// ----------------------- Setup RF Matrix ----------------------
627	MHMatrix train("Train");
628	train.AddColumns(fRules);
629	if (fEnableWeights[kTrainOn] \|\| fEnableWeights[kTrainOff])
630	train.AddColumn("MWeight.fVal");
631	train.AddColumn("MHadronness.fVal");
632
633	// ----------------------- Fill Matrix RF ----------------------
634
635	// Setup the hadronness container identifying gammas and off-data
636	// and setup a container for the weights
637	MParameterD had("MHadronness");
638	MParameterD wgt("MWeight");
639
640	// Add them to the parameter list
641	MParList plistx;
642	plistx.AddToList(this); // take care of fDisplay!
643	plistx.AddToList(&had);
644	plistx.AddToList(&wgt);
645
646	// Setup the tool class to fill the matrix
647	MTFillMatrix fill;
648	fill.SetLogStream(fLog);
649	fill.SetDisplay(fDisplay);
650	fill.AddPreCuts(fPreCuts);
651	fill.AddPreCuts(fTrainCuts);
652
653	// Set classifier for gammas
654	had.SetVal(0);
655	wgt.SetVal(1);
656
657	// Setup the tool class to read the gammas and read them
658	fill.SetName("FillGammas");
659	fill.SetDestMatrix1(&train, fNum[kTrainOn]);
660	fill.SetReader(&read1);
661	fill.AddPreTasks(fPreTasksSet[kTrainOn]);
662	fill.AddPreTasks(fPreTasks);
663	fill.AddPostTasks(fPostTasksSet[kTrainOn]);
664	fill.AddPostTasks(fPostTasks);
665	if (!fill.Process(plistx))
666	return kFALSE;
667
668	// Check the number or read events
669	const Int_t numgammastrn = train.GetNumRows();
670	if (numgammastrn==0)
671	{
672	*fLog << err << "ERROR - No gammas available for training... aborting." << endl;
673	return kFALSE;
674	}
675
676	// Remove possible post tasks
677	fill.ClearPreTasks();
678	fill.ClearPostTasks();
679
680	// Set classifier for background
681	had.SetVal(1);
682	wgt.SetVal(1);
683
684	// In case of wobble mode we have to do something special
685	MSrcPosRndm srcrndm;
686	srcrndm.SetDistOfSource(0.4);
687
688	MHillasCalc hcalc;
689	hcalc.SetFlags(MHillasCalc::kCalcHillasSrc);
690
691	if (fDataSetTrain.IsWobbleMode())
692	{
693	fPreTasksSet[kTrainOff].AddFirst(&hcalc);
694	fPreTasksSet[kTrainOff].AddFirst(&srcrndm);
695	}
696
697	// Setup the tool class to read the background and read them
698	fill.SetName("FillBackground");
699	fill.SetDestMatrix1(&train, fNum[kTrainOff]);
700	fill.SetReader(&read3);
701	fill.AddPreTasks(fPreTasksSet[kTrainOff]);
702	fill.AddPreTasks(fPreTasks);
703	fill.AddPostTasks(fPostTasksSet[kTrainOff]);
704	fill.AddPostTasks(fPostTasks);
705	if (!fill.Process(plistx))
706	return kFALSE;
707
708	// Check the number or read events
709	const Int_t numbackgrndtrn = train.GetNumRows()-numgammastrn;
710	if (numbackgrndtrn==0)
711	{
712	*fLog << err << "ERROR - No background available for training... aborting." << endl;
713	return kFALSE;
714	}
715
716	// ------------------------ Train RF --------------------------
717
718	MRanForestCalc rf;
719	rf.SetNumTrees(fNumTrees);
720	rf.SetNdSize(fNdSize);
721	rf.SetNumTry(fNumTry);
722	rf.SetNumObsoleteVariables(1);
723	rf.SetLastDataColumnHasWeights(fEnableWeights[kTrainOn] \|\| fEnableWeights[kTrainOff]);
724	rf.SetDebug(fDebug);
725	rf.SetDisplay(fDisplay);
726	rf.SetLogStream(fLog);
727	rf.SetFileName(out);
728	rf.SetNameOutput("MHadronness");
729
730	// Train the random forest either by classification or regression
731	if (fUseRegression)
732	{
733	if (!rf.TrainRegression(train)) // regression
734	return kFALSE;
735	}
736	else
737	{
738	if (!rf.TrainSingleRF(train)) // classification
739	return kFALSE;
740	}
741
742	// Output information about what was going on so far.
743	*fLog << all;
744	fLog->Separator("The forest was trained with...");
745
746	*fLog << "Training method:" << endl;
747	fLog << " " << (fUseRegression?"regression":"classification") << endl;
748	if (fEnableWeights[kTrainOn])
749	fLog << " weights for on-data" << endl;
750	if (fEnableWeights[kTrainOff])
751	fLog << " weights for off-data" << endl;
752	if (fDataSetTrain.IsWobbleMode())
753	fLog << " random source position in a distance of 0.4°" << endl;
754	*fLog << endl;
755	*fLog << "Events used for training:" << endl;
756	fLog << " Gammas: " << numgammastrn << endl;
757	fLog << " Background: " << numbackgrndtrn << endl;
758	*fLog << endl;
759	*fLog << "Gamma/Background ratio:" << endl;
760	fLog << " Requested: " << (float)fNum[kTrainOn]/fNum[kTrainOff] << endl;
761	fLog << " Result: " << (float)numgammastrn/numbackgrndtrn << endl;
762	*fLog << endl;
763	*fLog << "Run-Time: " << Form("%.1f", clock.RealTime()/60) << "min (CPU: ";
764	*fLog << Form("%.1f", clock.CpuTime()/60) << "min)" << endl;
765
766	// Chekc if testing is requested
767	if (!fDataSetTest.IsValid())
768	return kTRUE;
769
770	// --------------------- Display result ----------------------
771	fLog->Separator("Test");
772
773	clock.Continue();
774
775	// Setup parlist and tasklist for testing
776	MParList plist;
777	MTaskList tlist;
778	plist.AddToList(this); // Take care of display
779	plist.AddToList(&tlist);
780
781	MMcEvt mcevt;
782	plist.AddToList(&mcevt);
783
784	plist.AddToList(&wgt);
785
786	// ----- Setup histograms -----
787	MBinning binsy(50, 0 , 1, "BinningMH3Y", "lin");
788	MBinning binsx(40, 10, 100000, "BinningMH3X", "log");
789
790	plist.AddToList(&binsx);
791	plist.AddToList(&binsy);
792
793	MH3 h31("MHillas.fSize", "MHadronness.fVal");
794	MH3 h32("MHillas.fSize", "MHadronness.fVal");
795	MH3 h40("MMcEvt.fEnergy", "MHadronness.fVal");
796	h31.SetTitle("Background probability vs. Size:Size [phe]:Hadronness h");
797	h32.SetTitle("Background probability vs. Size:Size [phe]:Hadronness h");
798	h40.SetTitle("Background probability vs. Energy:Energy [GeV]:Hadronness h");
799
800	MHHadronness hist;
801
802	// ----- Setup tasks -----
803	MFillH fillh0(&hist, "", "FillHadronness");
804	MFillH fillh1(&h31);
805	MFillH fillh2(&h32);
806	MFillH fillh4(&h40);
807	fillh0.SetWeight("MWeight");
808	fillh1.SetWeight("MWeight");
809	fillh2.SetWeight("MWeight");
810	fillh4.SetWeight("MWeight");
811	fillh1.SetDrawOption("colz profy");
812	fillh2.SetDrawOption("colz profy");
813	fillh4.SetDrawOption("colz profy");
814	fillh1.SetNameTab("Background");
815	fillh2.SetNameTab("GammasH");
816	fillh4.SetNameTab("GammasE");
817	fillh0.SetBit(MFillH::kDoNotDisplay);
818
819	// ----- Setup filter -----
820	MFilterList precuts;
821	precuts.AddToList(fPreCuts);
822	precuts.AddToList(fTestCuts);
823
824	MContinue c0(&precuts);
825	c0.SetName("PreCuts");
826	c0.SetInverted();
827
828	MFEventSelector sel; // FIXME: USING IT (WITH PROB?) in READ will by much faster!!!
829	sel.SetNumSelectEvts(fNum[kTestOff]);
830
831	MContinue c1(&sel);
832	c1.SetInverted();
833
834	// ----- Setup tasklist -----
835	tlist.AddToList(&read2);
836	tlist.AddToList(&c1);
837	tlist.AddToList(fPreTasksSet[kTestOff]);
838	tlist.AddToList(fPreTasks);
839	tlist.AddToList(&c0);
840	tlist.AddToList(&rf);
841	tlist.AddToList(fPostTasksSet[kTestOff]);
842	tlist.AddToList(fPostTasks);
843	tlist.AddToList(&fillh0);
844	tlist.AddToList(&fillh1);
845
846	// Enable Acceleration
847	tlist.SetAccelerator(MTask::kAccDontReset\|MTask::kAccDontTime);
848
849	// ----- Run eventloop on background -----
850	MEvtLoop loop;
851	loop.SetDisplay(fDisplay);
852	loop.SetLogStream(fLog);
853	loop.SetParList(&plist);
854
855	wgt.SetVal(1);
856	if (!loop.Eventloop())
857	return kFALSE;
858	/*
859	if (!loop.GetDisplay())
860	{
861	gLog << warn << "Display closed by user... execution aborted." << endl << endl;
862	return kFALSE;
863	}
864	*/
865	// ----- Setup and run eventloop on gammas -----
866	sel.SetNumSelectEvts(fNum[kTestOn]);
867	fillh0.ResetBit(MFillH::kDoNotDisplay);
868
869	// Remove PreTasksOff and PostTasksOff from the list
870	tlist.RemoveFromList(fPreTasksSet[kTestOff]);
871	tlist.RemoveFromList(fPostTasksSet[kTestOff]);
872
873	// replace the reading task by a new one
874	tlist.Replace(&read4);
875
876	// Add the PreTasksOn directly after the reading task
877	tlist.AddToListAfter(fPreTasksSet[kTestOn], &c1);
878
879	// Add the PostTasksOn after rf
880	tlist.AddToListAfter(fPostTasksSet[kTestOn], &rf);
881
882	// Replace fillh1 by fillh2
883	tlist.Replace(&fillh2);
884
885	// Add fillh4 after the new fillh2
886	tlist.AddToListAfter(&fillh4, &fillh2);
887
888	// Enable Acceleration
889	tlist.SetAccelerator(MTask::kAccDontReset\|MTask::kAccDontTime);
890
891	wgt.SetVal(1);
892	if (!loop.Eventloop())
893	return kFALSE;
894
895	// Show what was going on in the testing
896	const Double_t numgammastst = h32.GetHist().GetEntries();
897	const Double_t numbackgrndtst = h31.GetHist().GetEntries();
898
899	*fLog << all;
900	fLog->Separator("The forest was tested with...");
901	*fLog << "Test method:" << endl;
902	fLog << " Random Forest: " << out << endl;
903	if (fEnableWeights[kTestOn])
904	fLog << " weights for on-data" << endl;
905	if (fEnableWeights[kTestOff])
906	fLog << " weights for off-data" << endl;
907	if (fDataSetTrain.IsWobbleMode())
908	fLog << " random source position in a distance of 0.4°" << endl;
909	*fLog << endl;
910	*fLog << "Events used for test:" << endl;
911	fLog << " Gammas: " << numgammastst << endl;
912	fLog << " Background: " << numbackgrndtst << endl;
913	*fLog << endl;
914	*fLog << "Gamma/Background ratio:" << endl;
915	fLog << " Requested: " << (float)fNum[kTestOn]/fNum[kTestOff] << endl;
916	fLog << " Result: " << (float)numgammastst/numbackgrndtst << endl;
917	*fLog << endl;
918
919	// Display the result plots
920	DisplayResult(h31, h32, ontime);
921
922	*fLog << "Total Run-Time: " << Form("%.1f", clock.RealTime()/60) << "min (CPU: ";
923	*fLog << Form("%.1f", clock.CpuTime()/60) << "min)" << endl;
924	fLog->Separator();
925
926	// Write the display
927	if (!WriteDisplay(out))
928	return kFALSE;
929
930	return kTRUE;
931	}
932

Note: See TracBrowser for help on using the repository browser.

Download in other formats: