Context Navigation

source: trunk/MagicSoft/Mars/mjobs/MJSpectrum.cc@ 6994

Visit:

Last change on this file since 6994 was 6994, checked in by tbretz, 20 years ago
* empty log message *
File size: 23.2 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, 4/2005 <mailto:tbretz@astro.uni-wuerzburg.de>
19	!
20	! Copyright: MAGIC Software Development, 2000-2005
21	!
22	!
23	\* ======================================================================== */
24
25	/////////////////////////////////////////////////////////////////////////////
26	//
27	// MJSpectrum
28	//
29	// Program to calculate spectrum
30	//
31	/////////////////////////////////////////////////////////////////////////////
32	#include "MJSpectrum.h"
33
34	// Root
35	#include <TF1.h>
36	#include <TH1.h>
37	#include <TH2.h>
38	#include <TFile.h>
39	#include <TChain.h>
40	#include <TCanvas.h>
41	#include <TObjArray.h>
42
43	// Environment
44	#include "MLog.h"
45	#include "MLogManip.h"
46
47	#include "MStatusArray.h"
48	#include "MStatusDisplay.h"
49
50	// Container
51	#include "MH3.h"
52	#include "MBinning.h"
53	#include "MDataSet.h"
54
55	// Spectrum
56	#include "../mhflux/MAlphaFitter.h"
57	#include "../mhflux/MHAlpha.h"
58	#include "../mhflux/MHCollectionArea.h"
59	#include "../mhflux/MHEnergyEst.h"
60
61	// Eventloop
62	#include "MEvtLoop.h"
63	#include "MTaskList.h"
64	#include "MParList.h"
65
66	// Tasks/Filter
67	#include "MReadMarsFile.h"
68	#include "MReadMarsFile.h"
69	#include "MFEventSelector2.h"
70	#include "MFDataMember.h"
71	#include "MEnergyEstimate.h"
72	#include "MTaskEnv.h"
73	#include "MFillH.h"
74	#include "MHillasCalc.h"
75	#include "MSrcPosCalc.h"
76	#include "MContinue.h"
77
78	ClassImp(MJSpectrum);
79
80	using namespace std;
81
82	MJSpectrum::MJSpectrum(const char name, const char title)
83	: fCut0(0),fCut1(0), fCut2(0), fCut3(0), fEstimateEnergy(0),
84	fRefill(kFALSE), fSimpleMode(kTRUE), fRawMc(kFALSE)
85	{
86	fName = name ? name : "MJSpectrum";
87	fTitle = title ? title : "Standard program to calculate spectrum";
88	}
89
90	MJSpectrum::~MJSpectrum()
91	{
92	if (fCut0)
93	delete fCut0;
94	if (fCut1)
95	delete fCut1;
96	if (fCut2)
97	delete fCut2;
98	if (fCut3)
99	delete fCut3;
100	if (fEstimateEnergy)
101	delete fEstimateEnergy;
102	}
103
104	// --------------------------------------------------------------------------
105	//
106	// Setup a task estimating the energy. The given task is cloned.
107	//
108	void MJSpectrum::SetEnergyEstimator(const MTask *task)
109	{
110	if (fEstimateEnergy)
111	delete fEstimateEnergy;
112	fEstimateEnergy = task ? (MTask*)task->Clone() : 0;
113	}
114
115	Bool_t MJSpectrum::ReadTask(MTask* &task, const char *name) const
116	{
117	if (task)
118	{
119	delete task;
120	task = 0;
121	}
122
123	task = (MTask*)gFile->Get(name);
124	if (!task)
125	{
126	*fLog << err << dbginf << "ERROR - " << name << " doen't exist in file!" << endl;
127	return kFALSE;
128	}
129	if (!task->InheritsFrom(MTask::Class()))
130	{
131	*fLog << err << dbginf << "ERROR - " << name << " read doesn't inherit from MTask!" << endl;
132	delete task;
133	return kFALSE;
134	}
135
136	task->SetName(name);
137
138	return kTRUE;
139	}
140
141	void MJSpectrum::PrintSetup(const MAlphaFitter &fit) const
142	{
143	fLog->Separator("Alpha Fitter");
144	*fLog << all;
145	fit.Print();
146
147	fLog->Separator("Used Cuts");
148	fCut0->Print();
149	fCut1->Print();
150	fCut2->Print();
151	fCut3->Print();
152
153	//gLog.Separator("Energy Estimator");
154	//fEstimateEnergy->Print();
155	}
156
157	Float_t MJSpectrum::ReadInput(MParList &plist, TH1D &h1, TH1D &h2)
158	{
159	*fLog << inf << "Reading from file: " << fPathIn << endl;
160
161	TFile file(fPathIn, "READ");
162	if (!file.IsOpen())
163	{
164	*fLog << err << dbginf << "ERROR - Could not open file " << fPathIn << endl;
165	return -1;
166	}
167
168	MStatusArray arr;
169	if (arr.Read()<=0)
170	{
171	*fLog << "MStatusDisplay not found in file... abort." << endl;
172	return -1;
173	}
174
175	TH1D vstime = (TH1D)arr.FindObjectInCanvas("Theta", "TH1D", "OnTime");
176	TH1D size = (TH1D)arr.FindObjectInCanvas("Excess", "TH1D", "Hist");
177	if (!vstime \|\| !size)
178	return -1;
179
180	vstime->Copy(h1);
181	size->Copy(h2);
182	h1.SetDirectory(0);
183	h2.SetDirectory(0);
184
185	if (fDisplay)
186	arr.DisplayIn(*fDisplay, "Hist");
187
188	if (!ReadTask(fCut0, "Cut0"))
189	return -1;
190	if (!ReadTask(fCut1, "Cut1"))
191	return -1;
192	if (!ReadTask(fCut2, "Cut2"))
193	return -1;
194	if (!ReadTask(fCut3, "Cut3"))
195	return -1;
196
197	TObjArray arrread;
198
199	TIter Next(plist);
200	TObject *o=0;
201	while ((o=Next()))
202	if (o->InheritsFrom(MBinning::Class()))
203	arrread.Add(o);
204
205	arrread.Add(plist.FindObject("MAlphaFitter"));
206
207	if (!ReadContainer(arrread))
208	return -1;
209
210	return vstime->Integral();
211	}
212
213	Bool_t MJSpectrum::ReadOrigMCDistribution(const MDataSet &set, TH1 &h) const
214	{
215	// Some debug output
216	fLog->Separator("Compiling original MC distribution");
217
218	*fLog << inf << "Please stand by, this may take a while..." << flush;
219
220	if (fDisplay)
221	fDisplay->SetStatusLine1("Compiling MC distribution...");
222
223	// Create chain
224	TChain chain("OriginalMC");
225	set.AddFilesOn(chain);
226
227	// Prepare histogram
228	h.Reset();
229
230	// Fill histogram from chain
231	h.SetDirectory(gROOT);
232	if (h.InheritsFrom(TH2::Class()))
233	{
234	h.SetNameTitle("ThetaEMC", "Event-Distribution vs Theta and Energy for MC (produced)");
235	h.SetXTitle("\\Theta [\\circ]");
236	h.SetYTitle("E [GeV]");
237	h.SetZTitle("Counts");
238	chain.Draw("MMcEvtBasic.fEnergy:MMcEvtBasic.fTelescopeTheta*TMath::RadToDeg()>>ThetaEMC", "", "goff");
239	}
240	else
241	{
242	h.SetNameTitle("ThetaMC", "Event-Distribution vs Theta for MC (produced)");
243	h.SetXTitle("\\Theta [\\circ]");
244	h.SetYTitle("Counts");
245	chain.Draw("MMcEvtBasic.fTelescopeTheta*TMath::RadToDeg()>>ThetaMC", "", "goff");
246	}
247	h.SetDirectory(0);
248
249	*fLog << "done." << endl;
250	if (fDisplay)
251	fDisplay->SetStatusLine2("done.");
252
253	if (h.GetEntries()>0)
254	return kTRUE;
255
256	*fLog << err << "ERROR - Histogram with original MC distribution empty..." << endl;
257
258	return h.GetEntries()>0;
259	}
260
261	Bool_t MJSpectrum::GetThetaDistribution(TH1D &temp1, TH1D &temp2) const
262	{
263	// Display some stuff
264	if (fDisplay)
265	{
266	TCanvas &c = fDisplay->AddTab("ZdDist");
267	c.Divide(2,2);
268
269	// On-Time vs. Theta
270	c.cd(1);
271	gPad->SetBorderMode(0);
272	temp1.DrawCopy();
273
274	// Number of MC events (produced) vs Theta
275	c.cd(2);
276	gPad->SetBorderMode(0);
277	temp2.SetName("NVsTheta");
278	temp2.DrawCopy();
279
280	c.cd(4);
281	gPad->SetBorderMode(0);
282
283	c.cd(3);
284	gPad->SetBorderMode(0);
285	}
286
287	// Calculate the Probability
288	temp1.Divide(&temp2);
289	temp1.Scale(1./temp1.GetMaximum());
290
291	// Some cosmetics: Name, Axis, etc.
292	temp1.SetName("ProbVsTheta");
293	temp1.SetTitle("Probability vs. Zenith Angle to choose MC events");
294	temp1.SetYTitle("Probability");
295	if (fDisplay)
296	temp1.DrawCopy();
297
298	return kTRUE;
299	}
300
301	void MJSpectrum::DisplayResult(const TH2D &h2) const
302	{
303	if (!fDisplay \|\| !fDisplay->CdCanvas("ZdDist"))
304	return;
305
306	TH1D &proj = *h2.ProjectionX();
307	proj.SetNameTitle("ThetaFinal", "Final Theta Distribution");
308	proj.SetXTitle("\\Theta [\\circ]");
309	proj.SetYTitle("Counts");
310	proj.SetLineColor(kBlue);
311	proj.SetDirectory(0);
312	proj.SetBit(kCanDelete);
313
314	TVirtualPad *pad = gPad;
315
316	pad->cd(4);
317	proj.DrawCopy();
318
319	pad->cd(1);
320	TH1D theta = (TH1D)gPad->FindObject("Theta");
321	if (theta)
322	{
323	proj.Scale(theta->GetMaximum()/proj.GetMaximum());
324	theta->SetMaximum(1.05*TMath::Max(theta->GetMaximum(), proj.GetMaximum()));
325	}
326	proj.Draw("same");
327	}
328
329	Bool_t MJSpectrum::Refill(MParList &plist, TH1D &h2) const
330	{
331	*fLog << endl;
332	fLog->Separator("Refill Excess");
333	*fLog << endl;
334
335	MTaskList tlist;
336	plist.AddToList(&tlist);
337
338	MReadTree read("Events");
339	read.DisableAutoScheme();
340	read.AddFile(fPathIn);
341
342	MEnergyEstimate est;
343	MTaskEnv taskenv1("EstimateEnergy");
344	taskenv1.SetDefault(fEstimateEnergy ? fEstimateEnergy : &est);
345
346	MFillH fill1("HistOff [MHAlpha]", "MHillasSrc", "FillHistOff");
347	MFillH fill2("HistOn [MHAlpha]", "MHillasSrc", "FillHistOn");
348
349	MFDataMember f0("DataType.fVal", '<', 0.5, "FilterOffData");
350	MFDataMember f1("DataType.fVal", '>', 0.5, "FilterOnData");
351
352	fill1.SetFilter(&f0);
353	fill2.SetFilter(&f1);
354
355	tlist.AddToList(&read);
356	tlist.AddToList(&taskenv1);
357	tlist.AddToList(&f0);
358	tlist.AddToList(&f1);
359	tlist.AddToList(&fill1);
360	tlist.AddToList(&fill2);
361
362	MEvtLoop loop(fName);
363	loop.SetParList(&plist);
364	loop.SetDisplay(fDisplay);
365	loop.SetLogStream(fLog);
366
367	if (!SetupEnv(loop))
368	return kFALSE;
369
370	if (!loop.Eventloop())
371	{
372	*fLog << err << GetDescriptor() << ": Refilling of data failed." << endl;
373	return kFALSE;
374	}
375
376	tlist.PrintStatistics();
377
378	if (!loop.GetDisplay())
379	{
380	*fLog << err << GetDescriptor() << ": Execution stopped by user." << endl;
381	return kFALSE;
382	}
383
384	const MHAlpha halpha = (MHAlpha )plist.FindObject("HistOn");
385	if (!halpha)
386	{
387	*fLog << err << GetDescriptor() << ": HistOn [MHAlpha] not found... abort." << endl;
388	return kFALSE;
389	}
390
391	halpha->GetHEnergy().Copy(h2);
392	h2.SetDirectory(0);
393
394	return kTRUE;
395	}
396
397	Bool_t MJSpectrum::IntermediateLoop(MParList &plist, MH3 &mh1, TH1D &temp1, const MDataSet &set) const
398	{
399	MTaskList tlist1;
400	plist.Replace(&tlist1);
401
402	MReadMarsFile readmc("OriginalMC");
403	//readmc.DisableAutoScheme();
404	set.AddFilesOn(readmc);
405	readmc.EnableBranch("MMcEvtBasic.fTelescopeTheta");
406	readmc.EnableBranch("MMcEvtBasic.fEnergy");
407
408	mh1.SetLogy();
409	mh1.SetLogz();
410	mh1.SetName("ThetaE");
411
412	MFillH fill0(&mh1);
413	//fill0.SetDrawOption("projx only");
414
415	MBinning bins2 = (MBinning)plist.FindObject("BinningEnergyEst");
416	MBinning bins3 = (MBinning)plist.FindObject("BinningTheta");
417	if (bins2 && bins3)
418	{
419	bins2->SetName("BinningThetaEY");
420	bins3->SetName("BinningThetaEX");
421	}
422	tlist1.AddToList(&readmc);
423
424	temp1.SetXTitle("MMcEvtBasic.fTelescopeTheta*kRad2Deg");
425	MH3 mh3mc(temp1);
426
427	MFEventSelector2 sel1(mh3mc);
428	sel1.SetHistIsProbability();
429
430	fill0.SetFilter(&sel1);
431
432	if (!fRawMc)
433	tlist1.AddToList(&sel1);
434	tlist1.AddToList(&fill0);
435
436	MEvtLoop loop1(fName);
437	loop1.SetParList(&plist);
438	loop1.SetLogStream(fLog);
439	loop1.SetDisplay(fDisplay);
440
441	if (!SetupEnv(loop1))
442	return kFALSE;
443
444	if (!loop1.Eventloop(fMaxEvents))
445	{
446	*fLog << err << GetDescriptor() << ": Processing of MC-data failed." << endl;
447	return kFALSE;
448	}
449
450	tlist1.PrintStatistics();
451
452	if (!loop1.GetDisplay())
453	{
454	*fLog << err << GetDescriptor() << ": Execution stopped by user." << endl;
455	return kFALSE;
456	}
457
458	if (bins2 && bins3)
459	{
460	bins2->SetName("BinningEnergyEst");
461	bins3->SetName("BinningTheta");
462	}
463	return kTRUE;
464	}
465
466	void MJSpectrum::DisplaySpectrum(MHCollectionArea &area, TH1D &excess, MHEnergyEst &hest, Double_t ontime) const
467	{
468	TH1D collarea(area.GetHEnergy());
469	TH1D spectrum(excess);
470	TH1D weights;
471	hest.GetWeights(weights);
472
473	cout << "Effective On time: " << ontime << "s" << endl;
474
475	spectrum.SetDirectory(NULL);
476	spectrum.SetBit(kCanDelete);
477	spectrum.Scale(1./ontime);
478	spectrum.Divide(&collarea);
479	spectrum.SetNameTitle("Preliminary", "N/sm^{2} versus Energy (before unfolding)");
480	spectrum.SetYTitle("N/sm^{2}");
481
482	TCanvas &c1 = fDisplay->AddTab("Spectrum");
483	c1.Divide(2,2);
484	c1.cd(1);
485	gPad->SetBorderMode(0);
486	gPad->SetLogx();
487	gPad->SetLogy();
488	gPad->SetGridx();
489	gPad->SetGridy();
490	collarea.DrawCopy();
491
492	c1.cd(2);
493	gPad->SetBorderMode(0);
494	gPad->SetLogx();
495	gPad->SetLogy();
496	gPad->SetGridx();
497	gPad->SetGridy();
498	spectrum.DrawCopy();
499
500	c1.cd(3);
501	gPad->SetBorderMode(0);
502	gPad->SetLogx();
503	gPad->SetLogy();
504	gPad->SetGridx();
505	gPad->SetGridy();
506	weights.DrawCopy();
507
508	//spectrum.Divide(&weights);
509	spectrum.Multiply(&weights);
510	spectrum.SetNameTitle("Flux", "N/TeVsm^{2} versus Energy (after unfolding)");
511
512	for (int i=0; i<excess.GetNbinsX(); i++)
513	{
514	spectrum.SetBinContent(i+1, spectrum.GetBinContent(i+1)/spectrum.GetBinWidth(i+1)*1000);
515	spectrum.SetBinError(i+1, spectrum.GetBinError(i+1)/ spectrum.GetBinWidth(i+1)*1000);
516	}
517
518	c1.cd(4);
519	gPad->SetBorderMode(0);
520	gPad->SetLogx();
521	gPad->SetLogy();
522	gPad->SetGridx();
523	gPad->SetGridy();
524	spectrum.SetXTitle("E [GeV]");
525	spectrum.SetYTitle("N/TeVsm^{2}");
526	spectrum.DrawCopy();
527
528	TF1 f("f", "[1]*(x/1e3)^[0]", 50, 3e4);
529	f.SetParameter(0, -2.87);
530	f.SetParameter(1, 1.9e-6);
531	f.SetLineColor(kGreen);
532	spectrum.Fit(&f, "NI", "", 55, 2e4);
533	f.DrawCopy("same");
534
535	/*
536	TString str;
537	str += "(1.68#pm0.15)10^{-7}";
538	str += "(\\frac{E}{TeV})^{-2.59#pm0.06}";
539	str += "\\frac{ph}{TeVm^{2}s}";
540
541	TLatex tex;
542	tex.DrawLatex(2e2, 7e-5, str);
543	*/
544	}
545
546	Bool_t MJSpectrum::PlotSame(MStatusArray &arr, MParList &plist, const char name, const char tab, const char *plot) const
547	{
548	cout << name << endl;
549	TString same(name);
550	same += "Same";
551
552	TH1 h1 = (TH1)arr.FindObjectInCanvas(name, "TH1F", tab);
553	TH1 h2 = (TH1)arr.FindObjectInCanvas(same, "TH1F", tab);
554	if (!h1 \|\| !h2)
555	return kFALSE;
556
557	TObject *obj = plist.FindObject(plot);
558	if (!obj)
559	{
560	*fLog << warn << plot << " not in parameter list... skipping." << endl;
561	return kFALSE;
562	}
563
564	TH1 h3 = (TH1)obj->FindObject(name);
565	if (!h3)
566	{
567	*fLog << warn << name << " not found in " << plot << "... skipping." << endl;
568	return kFALSE;
569	}
570
571
572	const MAlphaFitter fit = (MAlphaFitter)plist.FindObject("MAlphaFitter");
573	const Double_t scale = fit ? fit->GetScaleFactor() : 1;
574
575	gPad->SetBorderMode(0);
576	h2->SetLineColor(kBlack);
577	h3->SetLineColor(kBlue);
578	h2->Add(h1, -scale);
579
580	h2->Scale(1./h2->Integral());
581	h3->Scale(1./h3->Integral());
582
583	h2->SetMaximum(1.05*TMath::Max(h2->GetMaximum(), h3->GetMaximum()));
584
585	h2 = h2->DrawCopy();
586	h3 = h3->DrawCopy("same");
587
588	// Don't do this on the original object!
589	h2->SetStats(kFALSE);
590	h3->SetStats(kFALSE);
591
592	return kTRUE;
593	}
594
595	Bool_t MJSpectrum::DisplaySize(MParList &plist) const
596	{
597	*fLog << inf << "Reading from file: " << fPathIn << endl;
598
599	cout << "Opening..." << endl;
600	TFile file(fPathIn, "READ");
601	if (!file.IsOpen())
602	{
603	*fLog << err << dbginf << "ERROR - Could not open file " << fPathIn << endl;
604	return kFALSE;
605	}
606
607	cout << "Reading..." << endl;
608	file.cd();
609	MStatusArray arr;
610	if (arr.Read()<=0)
611	{
612	*fLog << "MStatusDisplay not found in file... abort." << endl;
613	return kFALSE;
614	}
615
616	cout << "Searching..." << endl;
617
618	TH1 excess = (TH1D)arr.FindObjectInCanvas("Excess", "TH1D", "Hist");
619	if (!excess)
620	return kFALSE;
621
622	cout << "Displaying..." << endl;
623
624	// ------------------- Plot excess versus size -------------------
625
626	TCanvas &c = fDisplay->AddTab("Excess");
627	c.Divide(3,2);
628	c.cd(1);
629	gPad->SetBorderMode(0);
630	gPad->SetLogx();
631	gPad->SetLogy();
632	gPad->SetGridx();
633	gPad->SetGridy();
634
635	excess->SetTitle("Number of excess events vs Size (data, mc/blue)");
636	excess->Scale(1./excess->Integral());
637	excess = excess->DrawCopy();
638	// Don't do this on the original object!
639	excess->SetStats(kFALSE);
640
641	cout << "ExcessSize..." << endl;
642
643	TObject *o=0;
644	if ((o=plist.FindObject("ExcessSize")))
645	{
646	TH1 histsel = (TH1F)o->FindObject("");
647	if (histsel)
648	{
649	histsel->Scale(1./histsel->Integral());
650	histsel->SetLineColor(kBlue);
651	histsel->SetBit(kCanDelete);
652	histsel = histsel->DrawCopy("same");
653	// Don't do this on the original object!
654	histsel->SetStats(kFALSE);
655	}
656	}
657
658	cout << "Dist..." << endl;
659
660	// -------------- Comparison of Image Parameters --------------
661	c.cd(2);
662	PlotSame(arr, plist, "Dist", "HilSrc", "MHHilSrcMCPost");
663
664	cout << "Length..." << endl;
665
666	c.cd(3);
667	PlotSame(arr, plist, "Length", "PostCut", "MHHillasMCPost");
668
669	cout << "M3l..." << endl;
670
671	c.cd(4);
672	PlotSame(arr, plist, "M3l", "HilExt", "MHHilExtMCPost");
673
674	cout << "Conc1..." << endl;
675
676	c.cd(5);
677	PlotSame(arr, plist, "Conc1", "NewPar", "MHNewParMCPost");
678
679	cout << "Width..." << endl;
680
681	c.cd(6);
682	PlotSame(arr, plist, "Width", "PostCut", "MHHillasMCPost");
683
684	return kTRUE;
685	}
686
687	Bool_t MJSpectrum::Process(const MDataSet &set)
688	{
689	if (!set.IsValid())
690	{
691	*fLog << err << "ERROR - DataSet invalid!" << endl;
692	return kFALSE;
693	}
694
695	CheckEnv();
696
697	// --------------------------------------------------------------------------------
698
699	*fLog << inf;
700	fLog->Separator(GetDescriptor());
701	*fLog << "Compile Monte Carlo Sample (data set " << set.GetName() << ")" << endl;
702	*fLog << endl;
703
704	MBinning bins1("BinningAlpha");
705	MBinning bins2("BinningEnergyEst");
706	MBinning bins3("BinningTheta");
707	MBinning bins4("BinningFalseSource");
708	MBinning bins5("BinningWidth");
709	MBinning bins6("BinningLength");
710	MBinning bins7("BinningDist");
711	MBinning bins8("BinningMaxDist");
712
713	MAlphaFitter fit;
714
715	MParList plist;
716	plist.AddToList(&bins1);
717	plist.AddToList(&bins2);
718	plist.AddToList(&bins3);
719	plist.AddToList(&bins4);
720	plist.AddToList(&bins5);
721	plist.AddToList(&bins6);
722	plist.AddToList(&bins7);
723	plist.AddToList(&bins8);
724	plist.AddToList(&fit);
725
726	TH1D temp1, size;
727	const Float_t ontime = ReadInput(plist, temp1, size);
728	if (ontime<0)
729	{
730	*fLog << err << GetDescriptor() << ": Could not determin effective on time..." << endl;
731	return kFALSE;
732	}
733
734	PrintSetup(fit);
735	bins3.SetEdges(temp1, 'x');
736
737	TH1D temp2(temp1);
738	if (!ReadOrigMCDistribution(set, temp2))
739	return kFALSE;
740
741	if (!GetThetaDistribution(temp1, temp2))
742	return kFALSE;
743
744	TH1D excess;
745	if (!Refill(plist, excess))
746	return kFALSE;
747
748	TH2D hist;
749	MH3 mh1("MMcEvtBasic.fTelescopeTheta*kRad2Deg", "MMcEvtBasic.fEnergy");
750	if (fSimpleMode)
751	{
752	hist.UseCurrentStyle();
753	MH::SetBinning(&hist, &bins3/temp1.GetXaxis()/, &bins2/excess.GetXaxis()/);
754	if (!ReadOrigMCDistribution(set, hist))
755	return kFALSE;
756
757	if (!fRawMc)
758	{
759	for (int y=0; y<hist.GetNbinsY(); y++)
760	for (int x=0; x<hist.GetNbinsX(); x++)
761	hist.SetBinContent(x, y, hist.GetBinContent(x, y)*temp1.GetBinContent(x));
762	}
763	}
764	else
765	if (!IntermediateLoop(plist, mh1, temp1, set))
766	return kFALSE;
767
768	DisplayResult(fSimpleMode ? hist : (TH2D&)mh1.GetHist());
769
770	// ------------------------- Final loop --------------------------
771
772	*fLog << endl;
773	fLog->Separator("Calculate efficiencies");
774	*fLog << endl;
775
776	MTaskList tlist2;
777	plist.Replace(&tlist2);
778
779	MReadMarsFile read("Events");
780	read.DisableAutoScheme();
781	set.AddFilesOn(read);
782
783	// Selector to get correct (final) theta-distribution
784	temp1.SetXTitle("MPointingPos.fZd");
785
786	MH3 mh3(temp1);
787
788	MFEventSelector2 sel2(mh3);
789	sel2.SetHistIsProbability();
790
791	MContinue contsel(&sel2);
792	contsel.SetInverted();
793
794	// Get correct source position
795	MSrcPosCalc calc;
796
797	// Calculate corresponding Hillas parameters
798	MHillasCalc hcalc1;
799	MHillasCalc hcalc2("MHillasCalcAnti");
800	hcalc1.SetFlags(MHillasCalc::kCalcHillasSrc);
801	hcalc2.SetFlags(MHillasCalc::kCalcHillasSrc);
802	hcalc2.SetNameHillasSrc("MHillasSrcAnti");
803	hcalc2.SetNameSrcPosCam("MSrcPosAnti");
804
805	// Fill collection area and energy estimator (unfolding)
806	// Make sure to use the same binning for MHCollectionArea and MHEnergyEst
807	MHCollectionArea area;
808	area.SetHistAll(fSimpleMode ? hist : (TH2D&)mh1.GetHist());
809	MHEnergyEst hest;
810
811	MFillH fill3(&area, "", "FillCollectionArea");
812	MFillH fill4(&hest, "", "FillEnergyEst");
813
814	MH3 hsize("MHillas.fSize");
815	//MH3 henergy("MEnergyEst.fVal");
816	hsize.SetName("ExcessSize");
817	//henergy.SetName("EnergyEst");
818	MBinning bins(size, "BinningExcessSize");
819	plist.AddToList(&hsize);
820	//plist.AddToList(&henergy);
821	plist.AddToList(&bins);
822
823	MFillH fill1a("MHHillasMCPre [MHHillas]", "MHillas", "FillHillasPre");
824	MFillH fill2a("MHHillasMCPost [MHHillas]", "MHillas", "FillHillasPost");
825	MFillH fill3a("MHVsSizeMCPost [MHVsSize]", "MHillasSrc", "FillVsSizePost");
826	MFillH fill4a("MHHilExtMCPost [MHHillasExt]", "MHillasSrc", "FillHilExtPost");
827	MFillH fill5a("MHHilSrcMCPost [MHHillasSrc]", "MHillasSrc", "FillHilSrcPost");
828	MFillH fill6a("MHImgParMCPost [MHImagePar]", "MImagePar", "FillImgParPost");
829	MFillH fill7a("MHNewParMCPost [MHNewImagePar]", "MNewImagePar", "FillNewParPost");
830	MFillH fill8a("ExcessSize [MH3]", "", "FillExcessSize");
831	//MFillH fill9a("EnergyEst [MH3]", "", "FillExcessEEst");
832	fill1a.SetNameTab("PreCut");
833	fill2a.SetNameTab("PostCut");
834	fill3a.SetNameTab("VsSize");
835	fill4a.SetNameTab("HilExt");
836	fill5a.SetNameTab("HilSrc");
837	fill6a.SetNameTab("ImgPar");
838	fill7a.SetNameTab("NewPar");
839	fill8a.SetBit(MFillH::kDoNotDisplay);
840	//fill9a.SetBit(MFillH::kDoNotDisplay);
841
842	MEnergyEstimate est;
843	MTaskEnv taskenv1("EstimateEnergy");
844	taskenv1.SetDefault(fEstimateEnergy ? fEstimateEnergy : &est);
845
846	tlist2.AddToList(&read);
847	if (!fRawMc)
848	tlist2.AddToList(&contsel);
849	tlist2.AddToList(&calc);
850	tlist2.AddToList(&hcalc1);
851	tlist2.AddToList(&hcalc2);
852	tlist2.AddToList(&fill1a);
853	tlist2.AddToList(fCut0);
854	tlist2.AddToList(fCut1);
855	tlist2.AddToList(fCut2);
856	tlist2.AddToList(fCut3);
857	tlist2.AddToList(&taskenv1);
858	tlist2.AddToList(&fill3);
859	tlist2.AddToList(&fill4);
860	tlist2.AddToList(&fill2a);
861	tlist2.AddToList(&fill3a);
862	tlist2.AddToList(&fill4a);
863	tlist2.AddToList(&fill5a);
864	tlist2.AddToList(&fill6a);
865	tlist2.AddToList(&fill7a);
866	tlist2.AddToList(&fill8a);
867	//tlist2.AddToList(&fill9a);
868
869	MEvtLoop loop2(fName);
870	loop2.SetParList(&plist);
871	loop2.SetDisplay(fDisplay);
872	loop2.SetLogStream(fLog);
873
874	if (!SetupEnv(loop2))
875	return kFALSE;
876
877	if (!loop2.Eventloop(fMaxEvents))
878	{
879	*fLog << err << GetDescriptor() << ": Processing of MC-data failed." << endl;
880	return kFALSE;
881	}
882
883	tlist2.PrintStatistics();
884
885	if (!loop2.GetDisplay())
886	{
887	*fLog << err << GetDescriptor() << ": Execution stopped by user." << endl;
888	return kFALSE;
889	}
890
891	// -------------------------- Spectrum ----------------------------
892
893	DisplaySpectrum(area, excess, hest, ontime);
894	DisplaySize(plist);
895
896	if (!fPathOut.IsNull())
897	fDisplay->SaveAsRoot(fPathOut);
898
899	return kTRUE;
900	}

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

Download in other formats: