Context Navigation

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

Visit:

Last change on this file since 8223 was 8223, checked in by tbretz, 18 years ago
* empty log message *
File size: 30.9 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(kGreen);
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 sig1 = MMath::SignificanceLiMa(s+b, b);
230	const Float_t sig2 = s<1 ? 0 : MMath::SignificanceLiMa(s+b, b)*TMath::Log10(s);
231
232	gr3.SetPoint(y, h.GetYaxis()->GetBinLowEdge(y+2), sig1);
233	gr4.SetPoint(y, h.GetYaxis()->GetBinLowEdge(y+2), sig2);
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 (!fDataSetTrain.IsValid())
592	{
593	*fLog << err << "ERROR - DataSet for training invalid!" << endl;
594	return kFALSE;
595	}
596	if (!fDataSetTest.IsValid())
597	{
598	*fLog << err << "ERROR - DataSet for testing invalid!" << endl;
599	return kFALSE;
600	}
601
602	if (fDataSetTrain.IsWobbleMode()!=fDataSetTest.IsWobbleMode())
603	{
604	*fLog << err << "ERROR - Train- and Test-DataSet have different observation modes!" << endl;
605	return kFALSE;
606	}
607
608	TStopwatch clock;
609	clock.Start();
610
611	// ----------------------- Auto Train? ----------------------
612
613	Float_t ontime = -1;
614	if (fAutoTrain)
615	{
616	fLog->Separator("Auto-Training -- Train-Data");
617	if (AutoTrain(fDataSetTrain, kTrainOn, kTrainOff, fFluxTrain)<0)
618	return kFALSE;
619	fLog->Separator("Auto-Training -- Test-Data");
620	ontime = AutoTrain(fDataSetTest, kTestOn, kTestOff, fFluxTest);
621	if (ontime<0)
622	return kFALSE;
623	}
624
625	// --------------------- Setup files --------------------
626	MReadReports read1;//("Events");
627	MReadMarsFile read2("Events");
628	MReadMarsFile read3("Events");
629	MReadReports read4;//("Events");
630	//read1.DisableAutoScheme();
631	read2.DisableAutoScheme();
632	read3.DisableAutoScheme();
633	//read4.DisableAutoScheme();
634
635	read1.AddTree("Events", "MTime.", MReadReports::kMaster);
636	read4.AddTree("Events", "MTime.", MReadReports::kMaster);
637	read1.AddTree("Drive", MReadReports::kRequired);
638	read4.AddTree("Drive", MReadReports::kRequired);
639	read1.AddTree("Starguider", MReadReports::kRequired);
640	read4.AddTree("Starguider", MReadReports::kRequired);
641
642	// Setup four reading tasks with the on- and off-data of the two datasets
643	if (!fDataSetTrain.AddFilesOn(read1))
644	return kFALSE;
645	const Bool_t setrc1 = fDataSetTrain.IsWobbleMode() ?
646	fDataSetTrain.AddFilesOn(read3) : fDataSetTrain.AddFilesOff(read3);
647	const Bool_t setrc2 = fDataSetTest.IsWobbleMode() ?
648	fDataSetTest.AddFilesOn(read2) : fDataSetTest.AddFilesOff(read2);
649	if (!setrc1 \|\| !setrc2)
650	return kFALSE;
651	if (!fDataSetTest.AddFilesOn(read4))
652	return kFALSE;
653
654	// ----------------------- Setup RF Matrix ----------------------
655	MHMatrix train("Train");
656	train.AddColumns(fRules);
657	if (fEnableWeights[kTrainOn] \|\| fEnableWeights[kTrainOff])
658	train.AddColumn("MWeight.fVal");
659	train.AddColumn("MHadronness.fVal");
660
661	// ----------------------- Fill Matrix RF ----------------------
662
663	// Setup the hadronness container identifying gammas and off-data
664	// and setup a container for the weights
665	MParameterD had("MHadronness");
666	MParameterD wgt("MWeight");
667
668	// Add them to the parameter list
669	MParList plistx;
670	plistx.AddToList(this); // take care of fDisplay!
671	plistx.AddToList(&had);
672	plistx.AddToList(&wgt);
673
674	// Setup the tool class to fill the matrix
675	MTFillMatrix fill;
676	fill.SetLogStream(fLog);
677	fill.SetDisplay(fDisplay);
678	fill.AddPreCuts(fPreCuts);
679	fill.AddPreCuts(fTrainCuts);
680
681	// Set classifier for gammas
682	had.SetVal(0);
683	wgt.SetVal(1);
684
685	// How to get source position from off- and on-data?
686	MPointingPos source("MSourcePos");
687	MObservatory obs;
688	MSrcPosCalc scalc;
689	MSrcPosCorrect scor;
690	MHillasCalc hcalcw;
691	MHillasCalc hcalcw2;
692	MPointingDevCalc devcalc;
693	scalc.SetMode(MSrcPosCalc::kDefault); // kWobble for off-source
694	hcalcw.SetFlags(MHillasCalc::kCalcHillasSrc);
695	hcalcw2.SetFlags(MHillasCalc::kCalcHillasSrc);
696	hcalcw2.SetNameHillasSrc("MHillasSrcAnti");
697	hcalcw2.SetNameSrcPosCam("MSrcPosAnti");
698	if (fDataSetTrain.IsWobbleMode())
699	{
700	// *******************************************************************
701	// Possible source position (eg. Wobble Mode)
702	if (fDataSetTrain.HasSource())
703	{
704	if (!fDataSetTrain.GetSourcePos(source))
705	return -1;
706	*fLog << all;
707	source.Print("RaDec");
708	}
709	else
710	*fLog << all << "No source position applied..." << endl;
711
712	// La Palma Magic1
713	plistx.AddToList(&obs);
714	plistx.AddToList(&source);
715
716	TList tlist2;
717	tlist2.Add(&scalc);
718	tlist2.Add(&scor);
719	tlist2.Add(&hcalcw);
720	tlist2.Add(&hcalcw2);
721
722	devcalc.SetStreamId("Starguider");
723	tlist2.Add(&devcalc);
724
725	fill.AddPreTasks(tlist2);
726	// *******************************************************************
727	}
728
729	// Setup the tool class to read the gammas and read them
730	fill.SetName("FillGammas");
731	fill.SetDestMatrix1(&train, fNum[kTrainOn]);
732	fill.SetReader(&read1);
733	fill.AddPreTasks(fPreTasksSet[kTrainOn]);
734	fill.AddPreTasks(fPreTasks);
735	fill.AddPostTasks(fPostTasksSet[kTrainOn]);
736	fill.AddPostTasks(fPostTasks);
737	if (!fill.Process(plistx))
738	return kFALSE;
739
740	// Check the number or read events
741	const Int_t numgammastrn = train.GetNumRows();
742	if (numgammastrn==0)
743	{
744	*fLog << err << "ERROR - No gammas available for training... aborting." << endl;
745	return kFALSE;
746	}
747
748	// Remove possible post tasks
749	fill.ClearPreTasks();
750	fill.ClearPostTasks();
751
752	// Set classifier for background
753	had.SetVal(1);
754	wgt.SetVal(1);
755
756	// In case of wobble mode we have to do something special
757	MSrcPosRndm srcrndm;
758	srcrndm.SetDistOfSource(0.4);
759
760	MHillasCalc hcalc;
761	MHillasCalc hcalc2;
762	hcalc.SetFlags(MHillasCalc::kCalcHillasSrc);
763	hcalc2.SetFlags(MHillasCalc::kCalcHillasSrc);
764	hcalc2.SetNameHillasSrc("MHillasSrcAnti");
765	hcalc2.SetNameSrcPosCam("MSrcPosAnti");
766
767	if (fDataSetTrain.IsWobbleMode())
768	{
769	scalc.SetMode(MSrcPosCalc::kWobble); // kWobble for off-source
770	fPreTasksSet[kTrainOff].AddFirst(&hcalc2);
771	fPreTasksSet[kTrainOff].AddFirst(&hcalc);
772	//fPreTasksSet[kTrainOff].AddFirst(&srcrndm);
773	fPreTasksSet[kTrainOff].AddFirst(&scor);
774	fPreTasksSet[kTrainOff].AddFirst(&scalc);
775	}
776
777	// Setup the tool class to read the background and read them
778	fill.SetName("FillBackground");
779	fill.SetDestMatrix1(&train, fNum[kTrainOff]);
780	fill.SetReader(&read3);
781	fill.AddPreTasks(fPreTasksSet[kTrainOff]);
782	fill.AddPreTasks(fPreTasks);
783	fill.AddPostTasks(fPostTasksSet[kTrainOff]);
784	fill.AddPostTasks(fPostTasks);
785	if (!fill.Process(plistx))
786	return kFALSE;
787
788	// Check the number or read events
789	const Int_t numbackgrndtrn = train.GetNumRows()-numgammastrn;
790	if (numbackgrndtrn==0)
791	{
792	*fLog << err << "ERROR - No background available for training... aborting." << endl;
793	return kFALSE;
794	}
795
796	// ------------------------ Train RF --------------------------
797
798	MRanForestCalc rf;
799	rf.SetNumTrees(fNumTrees);
800	rf.SetNdSize(fNdSize);
801	rf.SetNumTry(fNumTry);
802	rf.SetNumObsoleteVariables(1);
803	rf.SetLastDataColumnHasWeights(fEnableWeights[kTrainOn] \|\| fEnableWeights[kTrainOff]);
804	rf.SetDebug(fDebug);
805	rf.SetDisplay(fDisplay);
806	rf.SetLogStream(fLog);
807	rf.SetFileName(out);
808	rf.SetNameOutput("MHadronness");
809
810	// Train the random forest either by classification or regression
811	if (fUseRegression)
812	{
813	if (!rf.TrainRegression(train)) // regression
814	return kFALSE;
815	}
816	else
817	{
818	if (!rf.TrainSingleRF(train)) // classification
819	return kFALSE;
820	}
821
822	// Output information about what was going on so far.
823	*fLog << all;
824	fLog->Separator("The forest was trained with...");
825
826	*fLog << "Training method:" << endl;
827	fLog << " " << (fUseRegression?"regression":"classification") << endl;
828	if (fEnableWeights[kTrainOn])
829	fLog << " weights for on-data" << endl;
830	if (fEnableWeights[kTrainOff])
831	fLog << " weights for off-data" << endl;
832	if (fDataSetTrain.IsWobbleMode())
833	fLog << " random source position in a distance of 0.4°" << endl;
834	*fLog << endl;
835	*fLog << "Events used for training:" << endl;
836	fLog << " Gammas: " << numgammastrn << endl;
837	fLog << " Background: " << numbackgrndtrn << endl;
838	*fLog << endl;
839	*fLog << "Gamma/Background ratio:" << endl;
840	fLog << " Requested: " << (float)fNum[kTrainOn]/fNum[kTrainOff] << endl;
841	fLog << " Result: " << (float)numgammastrn/numbackgrndtrn << endl;
842	*fLog << endl;
843	*fLog << "Run-Time: " << Form("%.1f", clock.RealTime()/60) << "min (CPU: ";
844	*fLog << Form("%.1f", clock.CpuTime()/60) << "min)" << endl;
845	*fLog << endl;
846	*fLog << "Output file name: " << out << endl;
847
848	// Chekc if testing is requested
849	if (!fDataSetTest.IsValid())
850	return kTRUE;
851
852	// --------------------- Display result ----------------------
853	fLog->Separator("Test");
854
855	clock.Continue();
856
857	// Setup parlist and tasklist for testing
858	MParList plist;
859	MTaskList tlist;
860	plist.AddToList(this); // Take care of display
861	plist.AddToList(&tlist);
862
863	MMcEvt mcevt;
864	plist.AddToList(&mcevt);
865
866	plist.AddToList(&wgt);
867
868	// ----- Setup histograms -----
869	MBinning binsy(50, 0 , 1, "BinningMH3Y", "lin");
870	MBinning binsx(40, 10, 100000, "BinningMH3X", "log");
871
872	plist.AddToList(&binsx);
873	plist.AddToList(&binsy);
874
875	MH3 h31("MHillas.fSize", "MHadronness.fVal");
876	MH3 h32("MHillas.fSize", "MHadronness.fVal");
877	MH3 h40("MMcEvt.fEnergy", "MHadronness.fVal");
878	h31.SetTitle("Background probability vs. Size:Size [phe]:Hadronness h");
879	h32.SetTitle("Background probability vs. Size:Size [phe]:Hadronness h");
880	h40.SetTitle("Background probability vs. Energy:Energy [GeV]:Hadronness h");
881
882	MHHadronness hist;
883
884	// ----- Setup tasks -----
885	MFillH fillh0(&hist, "", "FillHadronness");
886	MFillH fillh1(&h31);
887	MFillH fillh2(&h32);
888	MFillH fillh4(&h40);
889	fillh0.SetWeight("MWeight");
890	fillh1.SetWeight("MWeight");
891	fillh2.SetWeight("MWeight");
892	fillh4.SetWeight("MWeight");
893	fillh1.SetDrawOption("colz profy");
894	fillh2.SetDrawOption("colz profy");
895	fillh4.SetDrawOption("colz profy");
896	fillh1.SetNameTab("Background");
897	fillh2.SetNameTab("GammasH");
898	fillh4.SetNameTab("GammasE");
899	fillh0.SetBit(MFillH::kDoNotDisplay);
900
901	// ----- Setup filter -----
902	MFilterList precuts;
903	precuts.AddToList(fPreCuts);
904	precuts.AddToList(fTestCuts);
905
906	MContinue c0(&precuts);
907	c0.SetName("PreCuts");
908	c0.SetInverted();
909
910	MFEventSelector sel; // FIXME: USING IT (WITH PROB?) in READ will by much faster!!!
911	sel.SetNumSelectEvts(fNum[kTestOff]);
912
913	MContinue c1(&sel);
914	c1.SetInverted();
915
916	// ----- Setup tasklist -----
917	tlist.AddToList(&read2);
918	if (fDataSetTest.IsWobbleMode())
919	{
920	tlist.AddToList(&srcrndm);
921	tlist.AddToList(&hcalc);
922	}
923	tlist.AddToList(&c1);
924	tlist.AddToList(fPreTasksSet[kTestOff]);
925	tlist.AddToList(fPreTasks);
926	tlist.AddToList(&c0);
927	tlist.AddToList(&rf);
928	tlist.AddToList(fPostTasksSet[kTestOff]);
929	tlist.AddToList(fPostTasks);
930	tlist.AddToList(&fillh0);
931	tlist.AddToList(&fillh1);
932
933	// Enable Acceleration
934	tlist.SetAccelerator(MTask::kAccDontReset\|MTask::kAccDontTime);
935
936	// ----- Run eventloop on background -----
937	MEvtLoop loop;
938	loop.SetDisplay(fDisplay);
939	loop.SetLogStream(fLog);
940	loop.SetParList(&plist);
941
942	wgt.SetVal(1);
943	if (!loop.Eventloop())
944	return kFALSE;
945	/*
946	if (!loop.GetDisplay())
947	{
948	gLog << warn << "Display closed by user... execution aborted." << endl << endl;
949	return kFALSE;
950	}
951	*/
952	// ----- Setup and run eventloop on gammas -----
953	sel.SetNumSelectEvts(fNum[kTestOn]);
954	fillh0.ResetBit(MFillH::kDoNotDisplay);
955
956	// Remove PreTasksOff and PostTasksOff from the list
957	tlist.RemoveFromList(fPreTasksSet[kTestOff]);
958	tlist.RemoveFromList(fPostTasksSet[kTestOff]);
959
960	// replace the reading task by a new one
961	tlist.Replace(&read4);
962
963	if (fDataSetTest.IsWobbleMode())
964	{
965	// *******************************************************************
966	// Possible source position (eg. Wobble Mode)
967	if (fDataSetTest.HasSource())
968	{
969	if (!fDataSetTest.GetSourcePos(source))
970	return -1;
971	*fLog << all;
972	source.Print("RaDec");
973	}
974	else
975	*fLog << all << "No source position applied..." << endl;
976
977	// La Palma Magic1
978	plist.AddToList(&obs);
979	plist.AddToList(&source);
980
981	// How to get source position from off- and on-data?
982	tlist.AddToListAfter(&scalc, &read4);
983	tlist.AddToListAfter(&scor, &scalc);
984	tlist.AddToListAfter(&hcalcw, &scor);
985
986	tlist.AddToList(&devcalc, "Starguider");
987	// *******************************************************************
988	}
989
990	// Add the PreTasksOn directly after the reading task
991	tlist.AddToListAfter(fPreTasksSet[kTestOn], &c1);
992
993	// Add the PostTasksOn after rf
994	tlist.AddToListAfter(fPostTasksSet[kTestOn], &rf);
995
996	// Replace fillh1 by fillh2
997	tlist.Replace(&fillh2);
998
999	// Add fillh4 after the new fillh2
1000	tlist.AddToListAfter(&fillh4, &fillh2);
1001
1002	// Enable Acceleration
1003	tlist.SetAccelerator(MTask::kAccDontReset\|MTask::kAccDontTime);
1004
1005	wgt.SetVal(1);
1006	if (!loop.Eventloop())
1007	return kFALSE;
1008
1009	// Show what was going on in the testing
1010	const Double_t numgammastst = h32.GetHist().GetEntries();
1011	const Double_t numbackgrndtst = h31.GetHist().GetEntries();
1012
1013	*fLog << all;
1014	fLog->Separator("The forest was tested with...");
1015	*fLog << "Test method:" << endl;
1016	fLog << " Random Forest: " << out << endl;
1017	if (fEnableWeights[kTestOn])
1018	fLog << " weights for on-data" << endl;
1019	if (fEnableWeights[kTestOff])
1020	fLog << " weights for off-data" << endl;
1021	if (fDataSetTrain.IsWobbleMode())
1022	fLog << " random source position in a distance of 0.4°" << endl;
1023	*fLog << endl;
1024	*fLog << "Events used for test:" << endl;
1025	fLog << " Gammas: " << numgammastst << endl;
1026	fLog << " Background: " << numbackgrndtst << endl;
1027	*fLog << endl;
1028	*fLog << "Gamma/Background ratio:" << endl;
1029	fLog << " Requested: " << (float)fNum[kTestOn]/fNum[kTestOff] << endl;
1030	fLog << " Result: " << (float)numgammastst/numbackgrndtst << endl;
1031	*fLog << endl;
1032
1033	// Display the result plots
1034	DisplayResult(h31, h32, ontime);
1035
1036	*fLog << "Total Run-Time: " << Form("%.1f", clock.RealTime()/60) << "min (CPU: ";
1037	*fLog << Form("%.1f", clock.CpuTime()/60) << "min)" << endl;
1038	fLog->Separator();
1039
1040	// Write the display
1041	if (!WriteDisplay(out))
1042	return kFALSE;
1043
1044	return kTRUE;
1045	}
1046

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

Download in other formats: