Context Navigation

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

Visit:

Last change on this file since 7023 was 7005, checked in by tbretz, 21 years ago
* empty log message *
File size: 23.5 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	// Try to find the class used to determin the signal!
332	TString cls("MHAlpha");
333	if (fDisplay)
334	{
335	TCanvas *c = fDisplay->GetCanvas("Hist");
336	if (c)
337	{
338	TIter Next(c->GetListOfPrimitives());
339	TObject *obj=0;
340	while ((obj=Next()))
341	if (obj->InheritsFrom(MHAlpha::Class()))
342	break;
343	if (obj)
344	cls = obj->ClassName();
345	}
346	}
347
348	cout << "FOUND: "<< cls << endl;
349
350	// Now fill the histogram
351	*fLog << endl;
352	fLog->Separator("Refill Excess");
353	*fLog << endl;
354
355	MTaskList tlist;
356	plist.AddToList(&tlist);
357
358	MReadTree read("Events");
359	read.DisableAutoScheme();
360	read.AddFile(fPathIn);
361
362	MEnergyEstimate est;
363	MTaskEnv taskenv1("EstimateEnergy");
364	taskenv1.SetDefault(fEstimateEnergy ? fEstimateEnergy : &est);
365
366	MFillH fill1(Form("HistEOff [%s]", cls.Data()), "", "FillHistEOff");
367	MFillH fill2(Form("HistE [%s]", cls.Data()), "", "FillHistE");
368
369	MFDataMember f0("DataType.fVal", '<', 0.5, "FilterOffData");
370	MFDataMember f1("DataType.fVal", '>', 0.5, "FilterOnData");
371
372	fill1.SetFilter(&f0);
373	fill2.SetFilter(&f1);
374
375	tlist.AddToList(&read);
376	tlist.AddToList(&taskenv1);
377	tlist.AddToList(&f0);
378	tlist.AddToList(&f1);
379	tlist.AddToList(&fill1);
380	tlist.AddToList(&fill2);
381
382	MEvtLoop loop(fName);
383	loop.SetParList(&plist);
384	loop.SetDisplay(fDisplay);
385	loop.SetLogStream(fLog);
386
387	if (!SetupEnv(loop))
388	return kFALSE;
389
390	if (!loop.Eventloop())
391	{
392	*fLog << err << GetDescriptor() << ": Refilling of data failed." << endl;
393	return kFALSE;
394	}
395
396	tlist.PrintStatistics();
397
398	if (!loop.GetDisplay())
399	{
400	*fLog << err << GetDescriptor() << ": Execution stopped by user." << endl;
401	return kFALSE;
402	}
403
404	const MHAlpha halpha = (MHAlpha )plist.FindObject("HistE");
405	if (!halpha)
406	{
407	*fLog << err << GetDescriptor() << ": HistE [MHAlpha] not found... abort." << endl;
408	return kFALSE;
409	}
410
411	halpha->GetHEnergy().Copy(h2);
412	h2.SetDirectory(0);
413
414	return kTRUE;
415	}
416
417	Bool_t MJSpectrum::IntermediateLoop(MParList &plist, MH3 &mh1, TH1D &temp1, const MDataSet &set) const
418	{
419	MTaskList tlist1;
420	plist.Replace(&tlist1);
421
422	MReadMarsFile readmc("OriginalMC");
423	//readmc.DisableAutoScheme();
424	set.AddFilesOn(readmc);
425	readmc.EnableBranch("MMcEvtBasic.fTelescopeTheta");
426	readmc.EnableBranch("MMcEvtBasic.fEnergy");
427
428	mh1.SetLogy();
429	mh1.SetLogz();
430	mh1.SetName("ThetaE");
431
432	MFillH fill0(&mh1);
433	//fill0.SetDrawOption("projx only");
434
435	MBinning bins2 = (MBinning)plist.FindObject("BinningEnergyEst");
436	MBinning bins3 = (MBinning)plist.FindObject("BinningTheta");
437	if (bins2 && bins3)
438	{
439	bins2->SetName("BinningThetaEY");
440	bins3->SetName("BinningThetaEX");
441	}
442	tlist1.AddToList(&readmc);
443
444	temp1.SetXTitle("MMcEvtBasic.fTelescopeTheta*kRad2Deg");
445	MH3 mh3mc(temp1);
446
447	MFEventSelector2 sel1(mh3mc);
448	sel1.SetHistIsProbability();
449
450	fill0.SetFilter(&sel1);
451
452	if (!fRawMc)
453	tlist1.AddToList(&sel1);
454	tlist1.AddToList(&fill0);
455
456	MEvtLoop loop1(fName);
457	loop1.SetParList(&plist);
458	loop1.SetLogStream(fLog);
459	loop1.SetDisplay(fDisplay);
460
461	if (!SetupEnv(loop1))
462	return kFALSE;
463
464	if (!loop1.Eventloop(fMaxEvents))
465	{
466	*fLog << err << GetDescriptor() << ": Processing of MC-data failed." << endl;
467	return kFALSE;
468	}
469
470	tlist1.PrintStatistics();
471
472	if (!loop1.GetDisplay())
473	{
474	*fLog << err << GetDescriptor() << ": Execution stopped by user." << endl;
475	return kFALSE;
476	}
477
478	if (bins2 && bins3)
479	{
480	bins2->SetName("BinningEnergyEst");
481	bins3->SetName("BinningTheta");
482	}
483	return kTRUE;
484	}
485
486	void MJSpectrum::DisplaySpectrum(MHCollectionArea &area, TH1D &excess, MHEnergyEst &hest, Double_t ontime) const
487	{
488	TH1D collarea(area.GetHEnergy());
489	TH1D spectrum(excess);
490	TH1D weights;
491	hest.GetWeights(weights);
492
493	cout << "Effective On time: " << ontime << "s" << endl;
494
495	spectrum.SetDirectory(NULL);
496	spectrum.SetBit(kCanDelete);
497	spectrum.Scale(1./ontime);
498	spectrum.Divide(&collarea);
499	spectrum.SetNameTitle("Preliminary", "N/sm^{2} versus Energy (before unfolding)");
500	spectrum.SetYTitle("N/sm^{2}");
501
502	TCanvas &c1 = fDisplay->AddTab("Spectrum");
503	c1.Divide(2,2);
504	c1.cd(1);
505	gPad->SetBorderMode(0);
506	gPad->SetLogx();
507	gPad->SetLogy();
508	gPad->SetGridx();
509	gPad->SetGridy();
510	collarea.DrawCopy();
511
512	c1.cd(2);
513	gPad->SetBorderMode(0);
514	gPad->SetLogx();
515	gPad->SetLogy();
516	gPad->SetGridx();
517	gPad->SetGridy();
518	spectrum.DrawCopy();
519
520	c1.cd(3);
521	gPad->SetBorderMode(0);
522	gPad->SetLogx();
523	gPad->SetLogy();
524	gPad->SetGridx();
525	gPad->SetGridy();
526	weights.DrawCopy();
527
528	//spectrum.Divide(&weights);
529	spectrum.Multiply(&weights);
530	spectrum.SetNameTitle("Flux", "N/TeVsm^{2} versus Energy (after unfolding)");
531
532	for (int i=0; i<excess.GetNbinsX(); i++)
533	{
534	spectrum.SetBinContent(i+1, spectrum.GetBinContent(i+1)/spectrum.GetBinWidth(i+1)*1000);
535	spectrum.SetBinError(i+1, spectrum.GetBinError(i+1)/ spectrum.GetBinWidth(i+1)*1000);
536	}
537
538	c1.cd(4);
539	gPad->SetBorderMode(0);
540	gPad->SetLogx();
541	gPad->SetLogy();
542	gPad->SetGridx();
543	gPad->SetGridy();
544	spectrum.SetXTitle("E [GeV]");
545	spectrum.SetYTitle("N/TeVsm^{2}");
546	spectrum.DrawCopy();
547
548	TF1 f("f", "[1]*(x/1e3)^[0]", 50, 3e4);
549	f.SetParameter(0, -2.87);
550	f.SetParameter(1, 1.9e-6);
551	f.SetLineColor(kGreen);
552	spectrum.Fit(&f, "NI", "", 55, 2e4);
553	f.DrawCopy("same");
554
555	/*
556	TString str;
557	str += "(1.68#pm0.15)10^{-7}";
558	str += "(\\frac{E}{TeV})^{-2.59#pm0.06}";
559	str += "\\frac{ph}{TeVm^{2}s}";
560
561	TLatex tex;
562	tex.DrawLatex(2e2, 7e-5, str);
563	*/
564	}
565
566	Bool_t MJSpectrum::PlotSame(MStatusArray &arr, MParList &plist, const char name, const char tab, const char *plot) const
567	{
568	cout << name << endl;
569	TString same(name);
570	same += "Same";
571
572	TH1 h1 = (TH1)arr.FindObjectInCanvas(name, "TH1F", tab);
573	TH1 h2 = (TH1)arr.FindObjectInCanvas(same, "TH1F", tab);
574	if (!h1 \|\| !h2)
575	return kFALSE;
576
577	TObject *obj = plist.FindObject(plot);
578	if (!obj)
579	{
580	*fLog << warn << plot << " not in parameter list... skipping." << endl;
581	return kFALSE;
582	}
583
584	TH1 h3 = (TH1)obj->FindObject(name);
585	if (!h3)
586	{
587	*fLog << warn << name << " not found in " << plot << "... skipping." << endl;
588	return kFALSE;
589	}
590
591
592	const MAlphaFitter fit = (MAlphaFitter)plist.FindObject("MAlphaFitter");
593	const Double_t scale = fit ? fit->GetScaleFactor() : 1;
594
595	gPad->SetBorderMode(0);
596	h2->SetLineColor(kBlack);
597	h3->SetLineColor(kBlue);
598	h2->Add(h1, -scale);
599
600	h2->Scale(1./h2->Integral());
601	h3->Scale(1./h3->Integral());
602
603	h2->SetMaximum(1.05*TMath::Max(h2->GetMaximum(), h3->GetMaximum()));
604
605	h2 = h2->DrawCopy();
606	h3 = h3->DrawCopy("same");
607
608	// Don't do this on the original object!
609	h2->SetStats(kFALSE);
610	h3->SetStats(kFALSE);
611
612	return kTRUE;
613	}
614
615	Bool_t MJSpectrum::DisplaySize(MParList &plist) const
616	{
617	*fLog << inf << "Reading from file: " << fPathIn << endl;
618
619	TFile file(fPathIn, "READ");
620	if (!file.IsOpen())
621	{
622	*fLog << err << dbginf << "ERROR - Could not open file " << fPathIn << endl;
623	return kFALSE;
624	}
625
626	file.cd();
627	MStatusArray arr;
628	if (arr.Read()<=0)
629	{
630	*fLog << "MStatusDisplay not found in file... abort." << endl;
631	return kFALSE;
632	}
633
634	TH1 excess = (TH1D)arr.FindObjectInCanvas("Excess", "TH1D", "Hist");
635	if (!excess)
636	return kFALSE;
637
638	// ------------------- Plot excess versus size -------------------
639
640	TCanvas &c = fDisplay->AddTab("Excess");
641	c.Divide(3,2);
642	c.cd(1);
643	gPad->SetBorderMode(0);
644	gPad->SetLogx();
645	gPad->SetLogy();
646	gPad->SetGridx();
647	gPad->SetGridy();
648
649	excess->SetTitle("Number of excess events vs Size (data, mc/blue)");
650	excess->Scale(1./excess->Integral());
651	excess = excess->DrawCopy();
652	// Don't do this on the original object!
653	excess->SetStats(kFALSE);
654
655	TObject *o=0;
656	if ((o=plist.FindObject("ExcessSize")))
657	{
658	TH1 histsel = (TH1F)o->FindObject("");
659	if (histsel)
660	{
661	histsel->Scale(1./histsel->Integral());
662	histsel->SetLineColor(kBlue);
663	histsel->SetBit(kCanDelete);
664	histsel = histsel->DrawCopy("same");
665	// Don't do this on the original object!
666	histsel->SetStats(kFALSE);
667	}
668	}
669
670	// -------------- Comparison of Image Parameters --------------
671	c.cd(2);
672	PlotSame(arr, plist, "Dist", "HilSrc", "MHHilSrcMCPost");
673
674	c.cd(3);
675	PlotSame(arr, plist, "Length", "PostCut", "MHHillasMCPost");
676
677	c.cd(4);
678	PlotSame(arr, plist, "M3l", "HilExt", "MHHilExtMCPost");
679
680	c.cd(5);
681	PlotSame(arr, plist, "Conc1", "NewPar", "MHNewParMCPost");
682
683	c.cd(6);
684	PlotSame(arr, plist, "Width", "PostCut", "MHHillasMCPost");
685
686	return kTRUE;
687	}
688
689	Bool_t MJSpectrum::Process(const MDataSet &set)
690	{
691	if (!set.IsValid())
692	{
693	*fLog << err << "ERROR - DataSet invalid!" << endl;
694	return kFALSE;
695	}
696
697	CheckEnv();
698
699	// --------------------------------------------------------------------------------
700
701	*fLog << inf;
702	fLog->Separator(GetDescriptor());
703	*fLog << "Compile Monte Carlo Sample (data set " << set.GetName() << ")" << endl;
704	*fLog << endl;
705
706	MBinning bins1("BinningAlpha");
707	MBinning bins2("BinningEnergyEst");
708	MBinning bins3("BinningTheta");
709	MBinning bins4("BinningFalseSource");
710	MBinning bins5("BinningWidth");
711	MBinning bins6("BinningLength");
712	MBinning bins7("BinningDist");
713	MBinning bins8("BinningMaxDist");
714	MBinning bins9("BinningM3Long");
715	MBinning bins0("BinningConc1");
716
717	MAlphaFitter fit;
718
719	MParList plist;
720	plist.AddToList(&bins1);
721	plist.AddToList(&bins2);
722	plist.AddToList(&bins3);
723	plist.AddToList(&bins4);
724	plist.AddToList(&bins5);
725	plist.AddToList(&bins6);
726	plist.AddToList(&bins7);
727	plist.AddToList(&bins8);
728	plist.AddToList(&bins9);
729	plist.AddToList(&bins0);
730	plist.AddToList(&fit);
731
732	TH1D temp1, size;
733	const Float_t ontime = ReadInput(plist, temp1, size);
734	if (ontime<0)
735	{
736	*fLog << err << GetDescriptor() << ": Could not determin effective on time..." << endl;
737	return kFALSE;
738	}
739
740	PrintSetup(fit);
741	bins3.SetEdges(temp1, 'x');
742
743	TH1D temp2(temp1);
744	if (!ReadOrigMCDistribution(set, temp2))
745	return kFALSE;
746
747	if (!GetThetaDistribution(temp1, temp2))
748	return kFALSE;
749
750	TH1D excess;
751	if (!Refill(plist, excess))
752	return kFALSE;
753
754	TH2D hist;
755	MH3 mh1("MMcEvtBasic.fTelescopeTheta*kRad2Deg", "MMcEvtBasic.fEnergy");
756	if (fSimpleMode)
757	{
758	hist.UseCurrentStyle();
759	MH::SetBinning(&hist, &bins3/temp1.GetXaxis()/, &bins2/excess.GetXaxis()/);
760	if (!ReadOrigMCDistribution(set, hist))
761	return kFALSE;
762
763	if (!fRawMc)
764	{
765	for (int y=0; y<hist.GetNbinsY(); y++)
766	for (int x=0; x<hist.GetNbinsX(); x++)
767	hist.SetBinContent(x, y, hist.GetBinContent(x, y)*temp1.GetBinContent(x));
768	}
769	}
770	else
771	if (!IntermediateLoop(plist, mh1, temp1, set))
772	return kFALSE;
773
774	DisplayResult(fSimpleMode ? hist : (TH2D&)mh1.GetHist());
775
776	// ------------------------- Final loop --------------------------
777
778	*fLog << endl;
779	fLog->Separator("Calculate efficiencies");
780	*fLog << endl;
781
782	MTaskList tlist2;
783	plist.Replace(&tlist2);
784
785	MReadMarsFile read("Events");
786	read.DisableAutoScheme();
787	set.AddFilesOn(read);
788
789	// Selector to get correct (final) theta-distribution
790	temp1.SetXTitle("MPointingPos.fZd");
791
792	MH3 mh3(temp1);
793
794	MFEventSelector2 sel2(mh3);
795	sel2.SetHistIsProbability();
796
797	MContinue contsel(&sel2);
798	contsel.SetInverted();
799
800	// Get correct source position
801	MSrcPosCalc calc;
802
803	// Calculate corresponding Hillas parameters
804	MHillasCalc hcalc1;
805	MHillasCalc hcalc2("MHillasCalcAnti");
806	hcalc1.SetFlags(MHillasCalc::kCalcHillasSrc);
807	hcalc2.SetFlags(MHillasCalc::kCalcHillasSrc);
808	hcalc2.SetNameHillasSrc("MHillasSrcAnti");
809	hcalc2.SetNameSrcPosCam("MSrcPosAnti");
810
811	// Fill collection area and energy estimator (unfolding)
812	// Make sure to use the same binning for MHCollectionArea and MHEnergyEst
813	MHCollectionArea area;
814	area.SetHistAll(fSimpleMode ? hist : (TH2D&)mh1.GetHist());
815	MHEnergyEst hest;
816
817	MFillH fill3(&area, "", "FillCollectionArea");
818	MFillH fill4(&hest, "", "FillEnergyEst");
819
820	MH3 hsize("MHillas.fSize");
821	//MH3 henergy("MEnergyEst.fVal");
822	hsize.SetName("ExcessSize");
823	//henergy.SetName("EnergyEst");
824	MBinning bins(size, "BinningExcessSize");
825	plist.AddToList(&hsize);
826	//plist.AddToList(&henergy);
827	plist.AddToList(&bins);
828
829	MFillH fill1a("MHHillasMCPre [MHHillas]", "MHillas", "FillHillasPre");
830	MFillH fill2a("MHHillasMCPost [MHHillas]", "MHillas", "FillHillasPost");
831	MFillH fill3a("MHVsSizeMCPost [MHVsSize]", "MHillasSrc", "FillVsSizePost");
832	MFillH fill4a("MHHilExtMCPost [MHHillasExt]", "MHillasSrc", "FillHilExtPost");
833	MFillH fill5a("MHHilSrcMCPost [MHHillasSrc]", "MHillasSrc", "FillHilSrcPost");
834	MFillH fill6a("MHImgParMCPost [MHImagePar]", "MImagePar", "FillImgParPost");
835	MFillH fill7a("MHNewParMCPost [MHNewImagePar]", "MNewImagePar", "FillNewParPost");
836	MFillH fill8a("ExcessSize [MH3]", "", "FillExcessSize");
837	//MFillH fill9a("EnergyEst [MH3]", "", "FillExcessEEst");
838	fill1a.SetNameTab("PreCut");
839	fill2a.SetNameTab("PostCut");
840	fill3a.SetNameTab("VsSize");
841	fill4a.SetNameTab("HilExt");
842	fill5a.SetNameTab("HilSrc");
843	fill6a.SetNameTab("ImgPar");
844	fill7a.SetNameTab("NewPar");
845	fill8a.SetBit(MFillH::kDoNotDisplay);
846	//fill9a.SetBit(MFillH::kDoNotDisplay);
847
848	MEnergyEstimate est;
849	MTaskEnv taskenv1("EstimateEnergy");
850	taskenv1.SetDefault(fEstimateEnergy ? fEstimateEnergy : &est);
851
852	tlist2.AddToList(&read);
853	if (!fRawMc)
854	tlist2.AddToList(&contsel);
855	tlist2.AddToList(&calc);
856	tlist2.AddToList(&hcalc1);
857	tlist2.AddToList(&hcalc2);
858	tlist2.AddToList(&fill1a);
859	tlist2.AddToList(fCut0);
860	tlist2.AddToList(fCut1);
861	tlist2.AddToList(fCut2);
862	tlist2.AddToList(fCut3);
863	tlist2.AddToList(&taskenv1);
864	tlist2.AddToList(&fill3);
865	tlist2.AddToList(&fill4);
866	tlist2.AddToList(&fill2a);
867	tlist2.AddToList(&fill3a);
868	tlist2.AddToList(&fill4a);
869	tlist2.AddToList(&fill5a);
870	tlist2.AddToList(&fill6a);
871	tlist2.AddToList(&fill7a);
872	tlist2.AddToList(&fill8a);
873	//tlist2.AddToList(&fill9a);
874
875	MEvtLoop loop2(fName);
876	loop2.SetParList(&plist);
877	loop2.SetDisplay(fDisplay);
878	loop2.SetLogStream(fLog);
879
880	if (!SetupEnv(loop2))
881	return kFALSE;
882
883	if (!loop2.Eventloop(fMaxEvents))
884	{
885	*fLog << err << GetDescriptor() << ": Processing of MC-data failed." << endl;
886	return kFALSE;
887	}
888
889	tlist2.PrintStatistics();
890
891	if (!loop2.GetDisplay())
892	{
893	*fLog << err << GetDescriptor() << ": Execution stopped by user." << endl;
894	return kFALSE;
895	}
896
897	// -------------------------- Spectrum ----------------------------
898
899	DisplaySpectrum(area, excess, hest, ontime);
900	DisplaySize(plist);
901
902	if (!fPathOut.IsNull())
903	fDisplay->SaveAsRoot(fPathOut);
904
905	return kTRUE;
906	}

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

Download in other formats: