Context Navigation

source: tags/Mars-V0.9.3/mjobs/MJSpectrum.cc

Visit:

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

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

Download in other formats: