Context Navigation

source: trunk/Mars/mjtrain/MJTrainSeparation.cc@ 19662

Visit:

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

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

Download in other formats: