Context Navigation

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

Visit:

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

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

Download in other formats: