Context Navigation

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

Visit:

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

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

Download in other formats: