Context Navigation

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

Visit:

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

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

Download in other formats: