Context Navigation

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

Visit:

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

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

Download in other formats: