Context Navigation

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

Visit:

Last change on this file since 8170 was 8170, checked in by tbretz, 18 years ago
* empty log message *
File size: 27.1 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	if (ontime>0)
291	*fLog << all << "Observation Time: " << TMath::Nint(ontime/60) << "min" << endl;
292	*fLog << "Maximum Significance: " << Form("%.1f", sig1);
293	if (ontime>0)
294	*fLog << Form(" [%.1f/sqrt(h)]", sig1/TMath::Sqrt(ontime/3600));
295	*fLog << endl;
296
297	*fLog << "Significance: S=" << Form("%.1f", sig1) << " E=" << s1 << " B=" << b1 << " h<";
298	*fLog << Form("%.2f", res2.GetYaxis()->GetBinCenter(my1)) << " s>";
299	*fLog << Form("%3d", TMath::Nint(res2.GetXaxis()->GetBinCenter(mx1))) << endl;
300	fLog << "SignificanceLogE: S=" << Form("%.1f", sig2/TMath::Log10(s2)) << " E=" << s2 << " B=" << b2 << " h<";
301	*fLog << Form("%.2f", res2.GetYaxis()->GetBinCenter(my2)) << " s>";
302	*fLog << Form("%3d", TMath::Nint(res2.GetXaxis()->GetBinCenter(mx2))) << endl;
303	*fLog << endl;
304	}
305
306	/*
307	Bool_t MJSpectrum::InitWeighting(const MDataSet &set, MMcSpectrumWeight &w) const
308	{
309	fLog->Separator("Initialize energy weighting");
310
311	if (!CheckEnv(w))
312	{
313	*fLog << err << "ERROR - Reading resources for MMcSpectrumWeight failed." << endl;
314	return kFALSE;
315	}
316
317	TChain chain("RunHeaders");
318	set.AddFilesOn(chain);
319
320	MMcCorsikaRunHeader *h=0;
321	chain.SetBranchAddress("MMcCorsikaRunHeader.", &h);
322	chain.GetEntry(1);
323
324	if (!h)
325	{
326	*fLog << err << "ERROR - Couldn't read MMcCorsikaRunHeader from DataSet." << endl;
327	return kFALSE;
328	}
329
330	if (!w.Set(*h))
331	{
332	*fLog << err << "ERROR - Initializing MMcSpectrumWeight failed." << endl;
333	return kFALSE;
334	}
335
336	w.Print();
337	return kTRUE;
338	}
339
340	Bool_t MJSpectrum::ReadOrigMCDistribution(const MDataSet &set, TH1 &h, MMcSpectrumWeight &weight) const
341	{
342	// Some debug output
343	fLog->Separator("Compiling original MC distribution");
344
345	weight.SetNameMcEvt("MMcEvtBasic");
346	const TString w(weight.GetFormulaWeights());
347	weight.SetNameMcEvt();
348
349	*fLog << inf << "Using weights: " << w << endl;
350	*fLog << "Please stand by, this may take a while..." << flush;
351
352	if (fDisplay)
353	fDisplay->SetStatusLine1("Compiling MC distribution...");
354
355	// Create chain
356	TChain chain("OriginalMC");
357	set.AddFilesOn(chain);
358
359	// Prepare histogram
360	h.Reset();
361
362	// Fill histogram from chain
363	h.SetDirectory(gROOT);
364	if (h.InheritsFrom(TH2::Class()))
365	{
366	h.SetNameTitle("ThetaEMC", "Event-Distribution vs Theta and Energy for MC (produced)");
367	h.SetXTitle("\\Theta [\\circ]");
368	h.SetYTitle("E [GeV]");
369	h.SetZTitle("Counts");
370	chain.Draw("MMcEvtBasic.fEnergy:MMcEvtBasic.fTelescopeTheta*TMath::RadToDeg()>>ThetaEMC", w, "goff");
371	}
372	else
373	{
374	h.SetNameTitle("ThetaMC", "Event-Distribution vs Theta for MC (produced)");
375	h.SetXTitle("\\Theta [\\circ]");
376	h.SetYTitle("Counts");
377	chain.Draw("MMcEvtBasic.fTelescopeTheta*TMath::RadToDeg()>>ThetaMC", w, "goff");
378	}
379	h.SetDirectory(0);
380
381	*fLog << "done." << endl;
382	if (fDisplay)
383	fDisplay->SetStatusLine2("done.");
384
385	if (h.GetEntries()>0)
386	return kTRUE;
387
388	*fLog << err << "ERROR - Histogram with original MC distribution empty..." << endl;
389
390	return h.GetEntries()>0;
391	}
392	*/
393
394	Bool_t MJTrainSeparation::GetEventsProduced(MDataSet &set, Double_t &num, Double_t &min, Double_t &max) const
395	{
396	TChain chain("OriginalMC");
397	set.AddFilesOn(chain);
398
399	min = chain.GetMinimum("MMcEvtBasic.fEnergy");
400	max = chain.GetMaximum("MMcEvtBasic.fEnergy");
401
402	num = chain.GetEntries();
403
404	if (num<100)
405	*fLog << err << "ERROR - Less than 100 entries in OriginalMC-Tree of MC-Train-Data found." << endl;
406
407	return num>=100;
408	}
409
410	Double_t MJTrainSeparation::GetDataRate(MDataSet &set, Double_t &num) const
411	{
412	TChain chain1("Events");
413	set.AddFilesOff(chain1);
414
415	num = chain1.GetEntries();
416	if (num<100)
417	{
418	*fLog << err << "ERROR - Less than 100 entries in Events-Tree of Train-Data found." << endl;
419	return -1;
420	}
421
422	TChain chain("EffectiveOnTime");
423	set.AddFilesOff(chain);
424
425	chain.Draw("MEffectiveOnTime.fVal", "MEffectiveOnTime.fVal", "goff");
426
427	TH1 h = dynamic_cast<TH1>(gROOT->FindObject("htemp"));
428	if (!h)
429	{
430	*fLog << err << "ERROR - Weird things are happening (htemp not found)!" << endl;
431	return -1;
432	}
433
434	const Double_t ontime = h->Integral();
435	delete h;
436
437	if (ontime<1)
438	{
439	*fLog << err << "ERROR - Less than 1s of effective observation time found in Train-Data." << endl;
440	return -1;
441	}
442
443	return num/ontime;
444	}
445
446	Double_t MJTrainSeparation::GetNumMC(MDataSet &set) const
447	{
448	TChain chain1("Events");
449	set.AddFilesOn(chain1);
450
451	const Double_t num = chain1.GetEntries();
452	if (num<100)
453	{
454	*fLog << err << "ERROR - Less than 100 entries in Events-Tree of Train-Data found." << endl;
455	return -1;
456	}
457
458	return num;
459	}
460
461	Float_t MJTrainSeparation::AutoTrain(MDataSet &set, Type_t typon, Type_t typof, Float_t flux)
462	{
463	Double_t num, min, max;
464	if (!GetEventsProduced(set, num, min, max))
465	return -1;
466
467	*fLog << inf << "Using build-in radius of 300m to calculate collection area!" << endl;
468
469	// Target spectrum
470	TF1 flx("Flux", "[0]/1000*(x/1000)^(-2.6)", min, max);
471	flx.SetParameter(0, flux);
472
473	// Number n0 of events this spectrum would produce per s and m^2
474	const Double_t n0 = flx.Integral(min, max); //[#]
475
476	// Area produced in MC
477	const Double_t A = TMath::Pi()300300; //[m²]
478
479	// Rate R of events this spectrum would produce per s
480	const Double_t R = n0*A; //[Hz]
481
482	fLog << "Source Spectrum: " << flux << " (E/TeV)^(-2.6) * TeVm^2s" << endl;
483
484	*fLog << "Gamma rate from the source inside the MC production area: " << R << "Hz" << endl;
485
486	// Number N of events produced (in trainings sample)
487	const Double_t N = num; //[#]
488
489	*fLog << "Events produced by MC inside the production area: " << TMath::Nint(num) << endl;
490
491	// This correponds to an observation time T [s]
492	const Double_t T = N/R; //[s]
493
494	*fLog << "Total time produced by the Monte Carlo: " << T << "s" << endl;
495
496	// With an average data rate after star of
497	Double_t data=0;
498	const Double_t r = GetDataRate(set, data); //[Hz]
499	Double_t ontime = data/r;
500
501	*fLog << "Events measured per second effective on time: " << r << "Hz" << endl;
502	*fLog << "Total effective on time: " << ontime << "s" << endl;
503
504	const Double_t ratio = T/ontime;
505	*fLog << "Ratio of Monte Carlo to data observation time: " << ratio << endl;
506
507	// 3570.5/43440.2 = 0.082
508
509
510	// this yields a number of n events to be read for training
511	const Double_t n = r*T; //[#]
512
513	*fLog << "Events to be read from the data sample: " << TMath::Nint(n) << endl;
514	*fLog << "Events available in data sample: " << data << endl;
515
516	if (r<0)
517	return -1;
518
519	Double_t nummc = GetNumMC(set);
520
521	*fLog << "Events available in MC sample: " << nummc << endl;
522
523	// *fLog << "MC read probability: " << data/n << endl;
524
525	// more data requested than available => Scale down num MC events
526	Double_t on, off;
527	if (data<n)
528	{
529	on = TMath::Nint(nummc*data/n);
530	off = TMath::Nint(data);
531	*fLog << warn;
532	*fLog << "Not enough data events available... scaling MC to data by " << data/n << endl;
533	*fLog << inf;
534	}
535	else
536	{
537	on = TMath::Nint(nummc);
538	off = TMath::Nint(n);
539	}
540
541	if (fNum[typon]>0 && fNum[typon]<on)
542	{
543	fNum[typof] = TMath::Nint(off*fNum[typon]/on);
544	ontime *= fNum[typon]/on;
545	*fLog << warn << "Less MC events requested... scaling data to MC by " << fNum[typon]/on << endl;
546	}
547	else
548	{
549	fNum[typon] = TMath::Nint(on);
550	fNum[typof] = TMath::Nint(off);
551	}
552
553	*fLog << inf;
554	*fLog << "Target number of MC events: " << fNum[typon] << endl;
555	*fLog << "Target number of data events: " << fNum[typof] << endl;
556
557	/*
558	An event rate dependent selection?
559	----------------------------------
560	Total average data rate: R
561	Goal number of events: N
562	Number of data events: N0
563	Rate assigned to single evt: r
564
565	Selection probability: N/N0 * r/R
566
567	f := N/N0 * r
568
569	MF f("f * MEventRate.fRate < rand");
570	*/
571
572	return ontime;
573	}
574
575	Bool_t MJTrainSeparation::Train(const char *out)
576	{
577	if (!fDataSetTrain.IsValid())
578	{
579	*fLog << err << "ERROR - DataSet for training invalid!" << endl;
580	return kFALSE;
581	}
582	if (!fDataSetTest.IsValid())
583	{
584	*fLog << err << "ERROR - DataSet for testing invalid!" << endl;
585	return kFALSE;
586	}
587
588	if (fDataSetTrain.IsWobbleMode()!=fDataSetTest.IsWobbleMode())
589	{
590	*fLog << err << "ERROR - Train- and Test-DataSet have different observation modes!" << endl;
591	return kFALSE;
592	}
593
594	TStopwatch clock;
595	clock.Start();
596
597	// ----------------------- Auto Train? ----------------------
598
599	Float_t ontime = -1;
600	if (fAutoTrain)
601	{
602	fLog->Separator("Auto-Training -- Train-Data");
603	if (AutoTrain(fDataSetTrain, kTrainOn, kTrainOff, fFluxTrain)<0)
604	return kFALSE;
605	fLog->Separator("Auto-Training -- Test-Data");
606	ontime = AutoTrain(fDataSetTest, kTestOn, kTestOff, fFluxTest);
607	if (ontime<0)
608	return kFALSE;
609	}
610
611	// --------------------- Setup files --------------------
612	MReadMarsFile read1("Events");
613	MReadMarsFile read2("Events");
614	MReadMarsFile read3("Events");
615	MReadMarsFile read4("Events");
616	read1.DisableAutoScheme();
617	read2.DisableAutoScheme();
618	read3.DisableAutoScheme();
619	read4.DisableAutoScheme();
620
621	// Setup four reading tasks with the on- and off-data of the two datasets
622	fDataSetTrain.AddFilesOn(read1);
623	fDataSetTrain.AddFilesOff(read3);
624
625	fDataSetTest.AddFilesOff(read2);
626	fDataSetTest.AddFilesOn(read4);
627
628	// ----------------------- Setup RF Matrix ----------------------
629	MHMatrix train("Train");
630	train.AddColumns(fRules);
631	if (fEnableWeights[kTrainOn] \|\| fEnableWeights[kTrainOff])
632	train.AddColumn("MWeight.fVal");
633	train.AddColumn("MHadronness.fVal");
634
635	// ----------------------- Fill Matrix RF ----------------------
636
637	// Setup the hadronness container identifying gammas and off-data
638	// and setup a container for the weights
639	MParameterD had("MHadronness");
640	MParameterD wgt("MWeight");
641
642	// Add them to the parameter list
643	MParList plistx;
644	plistx.AddToList(this); // take care of fDisplay!
645	plistx.AddToList(&had);
646	plistx.AddToList(&wgt);
647
648	// Setup the tool class to fill the matrix
649	MTFillMatrix fill;
650	fill.SetLogStream(fLog);
651	fill.SetDisplay(fDisplay);
652	fill.AddPreCuts(fPreCuts);
653	fill.AddPreCuts(fTrainCuts);
654
655	// Set classifier for gammas
656	had.SetVal(0);
657	wgt.SetVal(1);
658
659	// Setup the tool class to read the gammas and read them
660	fill.SetName("FillGammas");
661	fill.SetDestMatrix1(&train, fNum[kTrainOn]);
662	fill.SetReader(&read1);
663	fill.AddPreTasks(fPreTasksSet[kTrainOn]);
664	fill.AddPreTasks(fPreTasks);
665	fill.AddPostTasks(fPostTasksSet[kTrainOn]);
666	fill.AddPostTasks(fPostTasks);
667	if (!fill.Process(plistx))
668	return kFALSE;
669
670	// Check the number or read events
671	const Int_t numgammastrn = train.GetNumRows();
672	if (numgammastrn==0)
673	{
674	*fLog << err << "ERROR - No gammas available for training... aborting." << endl;
675	return kFALSE;
676	}
677
678	// Remove possible post tasks
679	fill.ClearPreTasks();
680	fill.ClearPostTasks();
681
682	// Set classifier for background
683	had.SetVal(1);
684	wgt.SetVal(1);
685
686	// In case of wobble mode we have to do something special
687	MSrcPosRndm srcrndm;
688	srcrndm.SetDistOfSource(0.4);
689
690	MHillasCalc hcalc;
691	hcalc.SetFlags(MHillasCalc::kCalcHillasSrc);
692
693	if (fDataSetTrain.IsWobbleMode())
694	{
695	fPreTasksSet[kTrainOff].AddFirst(&hcalc);
696	fPreTasksSet[kTrainOff].AddFirst(&srcrndm);
697	}
698
699	// Setup the tool class to read the background and read them
700	fill.SetName("FillBackground");
701	fill.SetDestMatrix1(&train, fNum[kTrainOff]);
702	fill.SetReader(&read3);
703	fill.AddPreTasks(fPreTasksSet[kTrainOff]);
704	fill.AddPreTasks(fPreTasks);
705	fill.AddPostTasks(fPostTasksSet[kTrainOff]);
706	fill.AddPostTasks(fPostTasks);
707	if (!fill.Process(plistx))
708	return kFALSE;
709
710	// Check the number or read events
711	const Int_t numbackgrndtrn = train.GetNumRows()-numgammastrn;
712	if (numbackgrndtrn==0)
713	{
714	*fLog << err << "ERROR - No background available for training... aborting." << endl;
715	return kFALSE;
716	}
717
718	// ------------------------ Train RF --------------------------
719
720	MRanForestCalc rf;
721	rf.SetNumTrees(fNumTrees);
722	rf.SetNdSize(fNdSize);
723	rf.SetNumTry(fNumTry);
724	rf.SetNumObsoleteVariables(1);
725	rf.SetLastDataColumnHasWeights(fEnableWeights[kTrainOn] \|\| fEnableWeights[kTrainOff]);
726	rf.SetDebug(fDebug);
727	rf.SetDisplay(fDisplay);
728	rf.SetLogStream(fLog);
729	rf.SetFileName(out);
730	rf.SetNameOutput("MHadronness");
731
732	// Train the random forest either by classification or regression
733	if (fUseRegression)
734	{
735	if (!rf.TrainRegression(train)) // regression
736	return kFALSE;
737	}
738	else
739	{
740	if (!rf.TrainSingleRF(train)) // classification
741	return kFALSE;
742	}
743
744	// Output information about what was going on so far.
745	*fLog << all;
746	fLog->Separator("The forest was trained with...");
747
748	*fLog << "Training method:" << endl;
749	fLog << " " << (fUseRegression?"regression":"classification") << endl;
750	if (fEnableWeights[kTrainOn])
751	fLog << " weights for on-data" << endl;
752	if (fEnableWeights[kTrainOff])
753	fLog << " weights for off-data" << endl;
754	if (fDataSetTrain.IsWobbleMode())
755	fLog << " random source position in a distance of 0.4°" << endl;
756	*fLog << endl;
757	*fLog << "Events used for training:" << endl;
758	fLog << " Gammas: " << numgammastrn << endl;
759	fLog << " Background: " << numbackgrndtrn << endl;
760	*fLog << endl;
761	*fLog << "Gamma/Background ratio:" << endl;
762	fLog << " Requested: " << (float)fNum[kTrainOn]/fNum[kTrainOff] << endl;
763	fLog << " Result: " << (float)numgammastrn/numbackgrndtrn << endl;
764	*fLog << endl;
765	*fLog << "Run-Time: " << Form("%.1f", clock.RealTime()/60) << "min (CPU: ";
766	*fLog << Form("%.1f", clock.CpuTime()/60) << "min)" << endl;
767	*fLog << "Output file name: " << out << endl;
768
769	// Chekc if testing is requested
770	if (!fDataSetTest.IsValid())
771	return kTRUE;
772
773	// --------------------- Display result ----------------------
774	fLog->Separator("Test");
775
776	clock.Continue();
777
778	// Setup parlist and tasklist for testing
779	MParList plist;
780	MTaskList tlist;
781	plist.AddToList(this); // Take care of display
782	plist.AddToList(&tlist);
783
784	MMcEvt mcevt;
785	plist.AddToList(&mcevt);
786
787	plist.AddToList(&wgt);
788
789	// ----- Setup histograms -----
790	MBinning binsy(50, 0 , 1, "BinningMH3Y", "lin");
791	MBinning binsx(40, 10, 100000, "BinningMH3X", "log");
792
793	plist.AddToList(&binsx);
794	plist.AddToList(&binsy);
795
796	MH3 h31("MHillas.fSize", "MHadronness.fVal");
797	MH3 h32("MHillas.fSize", "MHadronness.fVal");
798	MH3 h40("MMcEvt.fEnergy", "MHadronness.fVal");
799	h31.SetTitle("Background probability vs. Size:Size [phe]:Hadronness h");
800	h32.SetTitle("Background probability vs. Size:Size [phe]:Hadronness h");
801	h40.SetTitle("Background probability vs. Energy:Energy [GeV]:Hadronness h");
802
803	MHHadronness hist;
804
805	// ----- Setup tasks -----
806	MFillH fillh0(&hist, "", "FillHadronness");
807	MFillH fillh1(&h31);
808	MFillH fillh2(&h32);
809	MFillH fillh4(&h40);
810	fillh0.SetWeight("MWeight");
811	fillh1.SetWeight("MWeight");
812	fillh2.SetWeight("MWeight");
813	fillh4.SetWeight("MWeight");
814	fillh1.SetDrawOption("colz profy");
815	fillh2.SetDrawOption("colz profy");
816	fillh4.SetDrawOption("colz profy");
817	fillh1.SetNameTab("Background");
818	fillh2.SetNameTab("GammasH");
819	fillh4.SetNameTab("GammasE");
820	fillh0.SetBit(MFillH::kDoNotDisplay);
821
822	// ----- Setup filter -----
823	MFilterList precuts;
824	precuts.AddToList(fPreCuts);
825	precuts.AddToList(fTestCuts);
826
827	MContinue c0(&precuts);
828	c0.SetName("PreCuts");
829	c0.SetInverted();
830
831	MFEventSelector sel; // FIXME: USING IT (WITH PROB?) in READ will by much faster!!!
832	sel.SetNumSelectEvts(fNum[kTestOff]);
833
834	MContinue c1(&sel);
835	c1.SetInverted();
836
837	// ----- Setup tasklist -----
838	tlist.AddToList(&read2);
839	tlist.AddToList(&c1);
840	tlist.AddToList(fPreTasksSet[kTestOff]);
841	tlist.AddToList(fPreTasks);
842	tlist.AddToList(&c0);
843	tlist.AddToList(&rf);
844	tlist.AddToList(fPostTasksSet[kTestOff]);
845	tlist.AddToList(fPostTasks);
846	tlist.AddToList(&fillh0);
847	tlist.AddToList(&fillh1);
848
849	// Enable Acceleration
850	tlist.SetAccelerator(MTask::kAccDontReset\|MTask::kAccDontTime);
851
852	// ----- Run eventloop on background -----
853	MEvtLoop loop;
854	loop.SetDisplay(fDisplay);
855	loop.SetLogStream(fLog);
856	loop.SetParList(&plist);
857
858	wgt.SetVal(1);
859	if (!loop.Eventloop())
860	return kFALSE;
861	/*
862	if (!loop.GetDisplay())
863	{
864	gLog << warn << "Display closed by user... execution aborted." << endl << endl;
865	return kFALSE;
866	}
867	*/
868	// ----- Setup and run eventloop on gammas -----
869	sel.SetNumSelectEvts(fNum[kTestOn]);
870	fillh0.ResetBit(MFillH::kDoNotDisplay);
871
872	// Remove PreTasksOff and PostTasksOff from the list
873	tlist.RemoveFromList(fPreTasksSet[kTestOff]);
874	tlist.RemoveFromList(fPostTasksSet[kTestOff]);
875
876	// replace the reading task by a new one
877	tlist.Replace(&read4);
878
879	// Add the PreTasksOn directly after the reading task
880	tlist.AddToListAfter(fPreTasksSet[kTestOn], &c1);
881
882	// Add the PostTasksOn after rf
883	tlist.AddToListAfter(fPostTasksSet[kTestOn], &rf);
884
885	// Replace fillh1 by fillh2
886	tlist.Replace(&fillh2);
887
888	// Add fillh4 after the new fillh2
889	tlist.AddToListAfter(&fillh4, &fillh2);
890
891	// Enable Acceleration
892	tlist.SetAccelerator(MTask::kAccDontReset\|MTask::kAccDontTime);
893
894	wgt.SetVal(1);
895	if (!loop.Eventloop())
896	return kFALSE;
897
898	// Show what was going on in the testing
899	const Double_t numgammastst = h32.GetHist().GetEntries();
900	const Double_t numbackgrndtst = h31.GetHist().GetEntries();
901
902	*fLog << all;
903	fLog->Separator("The forest was tested with...");
904	*fLog << "Test method:" << endl;
905	fLog << " Random Forest: " << out << endl;
906	if (fEnableWeights[kTestOn])
907	fLog << " weights for on-data" << endl;
908	if (fEnableWeights[kTestOff])
909	fLog << " weights for off-data" << endl;
910	if (fDataSetTrain.IsWobbleMode())
911	fLog << " random source position in a distance of 0.4°" << endl;
912	*fLog << endl;
913	*fLog << "Events used for test:" << endl;
914	fLog << " Gammas: " << numgammastst << endl;
915	fLog << " Background: " << numbackgrndtst << endl;
916	*fLog << endl;
917	*fLog << "Gamma/Background ratio:" << endl;
918	fLog << " Requested: " << (float)fNum[kTestOn]/fNum[kTestOff] << endl;
919	fLog << " Result: " << (float)numgammastst/numbackgrndtst << endl;
920	*fLog << endl;
921
922	// Display the result plots
923	DisplayResult(h31, h32, ontime);
924
925	*fLog << "Total Run-Time: " << Form("%.1f", clock.RealTime()/60) << "min (CPU: ";
926	*fLog << Form("%.1f", clock.CpuTime()/60) << "min)" << endl;
927	fLog->Separator();
928
929	// Write the display
930	if (!WriteDisplay(out))
931	return kFALSE;
932
933	return kTRUE;
934	}
935

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

Download in other formats: