Context Navigation

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

Visit:

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

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

Download in other formats: