Context Navigation

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

Visit:

Last change on this file since 7600 was 7522, checked in by tbretz, 19 years ago
* empty log message *
File size: 33.0 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	TH1D spec=(TH1D)spectrum.DrawCopy();
665	FitSpectrum(*spec);
666
667	c1.cd(3);
668	gPad->SetBorderMode(0);
669	gPad->SetLogx();
670	gPad->SetLogy();
671	gPad->SetGridx();
672	gPad->SetGridy();
673	weights.DrawCopy();
674
675	//spectrum.Multiply(&weights);
676	spectrum.SetNameTitle("Flux", "Spectrum (after unfolding)");
677	spectrum.SetBit(TH1::kNoStats);
678
679	for (int i=0; i<excess.GetNbinsX(); i++)
680	{
681	spectrum.SetBinContent(i+1, spectrum.GetBinContent(i+1)*weights.GetBinContent(i+1));
682	spectrum.SetBinError(i+1, spectrum.GetBinError(i+1) *weights.GetBinContent(i+1));
683
684	spectrum.SetBinContent(i+1, spectrum.GetBinContent(i+1)/spectrum.GetBinWidth(i+1)*1000);
685	spectrum.SetBinError(i+1, spectrum.GetBinError(i+1)/ spectrum.GetBinWidth(i+1)*1000);
686	}
687
688	c1.cd(4);
689	gPad->SetBorderMode(0);
690	gPad->SetLogx();
691	gPad->SetLogy();
692	gPad->SetGridx();
693	gPad->SetGridy();
694	spectrum.SetMinimum(1e-12);
695	spectrum.SetXTitle("E [GeV]");
696	spectrum.SetYTitle("dN/dE [N/TeVsm^{2}]");
697	spec = (TH1D*)spectrum.DrawCopy();
698
699	// Calculate Spectrum * E^2
700	for (int i=0; i<spectrum.GetNbinsX(); i++)
701	{
702	const Double_t e = TMath::Sqrt(spectrum.GetBinLowEdge(i+1)spectrum.GetBinLowEdge(i+2))1e-3;
703
704	spectrum.SetBinContent(i+1, spectrum.GetBinContent(i+1)ee);
705	spectrum.SetBinError( i+1, spectrum.GetBinError(i+1) ee);
706	}
707
708	// Display dN/dE*E^2 for conveinience
709	fDisplay->AddTab("Check");
710	gPad->SetLogx();
711	gPad->SetLogy();
712	gPad->SetGrid();
713
714	spectrum.SetName("FluxStd");
715	spectrum.SetMarkerStyle(kFullDotMedium);
716	spectrum.SetTitle("Differential flux times E^{2}");
717	spectrum.SetYTitle("E^{2}dN/dE [N TeV/sm^{2}]");
718	spectrum.SetDirectory(0);
719	spectrum.DrawCopy();
720
721	// Fit Spectrum
722	c1.cd(4);
723	return FitSpectrum(*spec);
724
725	/*
726	TF1 f("f", "[1]*(x/1e3)^[0]", 10, 3e4);
727	f.SetParameter(0, -2.87);
728	f.SetParameter(1, 1.9e-6);
729	f.SetLineColor(kGreen);
730	spectrum.Fit(&f, "NIM", "", 100, 5000);
731	f.DrawCopy("same");
732
733	const Double_t p0 = f.GetParameter(0);
734	const Double_t p1 = f.GetParameter(1);
735
736	const Double_t e0 = f.GetParError(0);
737	const Double_t e1 = f.GetParError(1);
738
739	const Int_t np = TMath::Nint(TMath::Floor(TMath::Log10(p1)));
740	const Double_t exp = TMath::Power(10, np);
741
742	TString str;
743	str += Form("(%.2f#pm%.2f)10^{%d}", p1/exp, e1/exp, np);
744	str += Form("(\\frac{E}{TeV})^{%.2f#pm%.2f}", p0, e0);
745	str += "\\frac{ph}{TeVm^{2}s}";
746
747	TLatex tex;
748	tex.SetTextSize(0.045);
749	tex.SetBit(TLatex::kTextNDC);
750	tex.DrawLatex(0.45, 0.935, str);
751
752	str = Form("\\chi^{2}/NDF=%.2f", f.GetChisquare()/f.GetNDF());
753	tex.DrawLatex(0.70, 0.83, str);
754
755	TArrayD res(2);
756	res[0] = f.GetParameter(0);
757	res[1] = f.GetParameter(1);
758
759	return res;
760	*/
761	/*
762	// Plot other spectra from Whipple
763	f.SetParameter(0, -2.45);
764	f.SetParameter(1, 3.3e-7);
765	f.SetRange(300, 8000);
766	f.SetLineColor(kBlack);
767	f.SetLineStyle(kDashed);
768	f.DrawCopy("same");
769
770	// Plot other spectra from Cangaroo
771	f.SetParameter(0, -2.53);
772	f.SetParameter(1, 2.0e-7);
773	f.SetRange(7000, 50000);
774	f.SetLineColor(kBlack);
775	f.SetLineStyle(kDashed);
776	f.DrawCopy("same");
777
778	// Plot other spectra from Robert
779	f.SetParameter(0, -2.59);
780	f.SetParameter(1, 2.58e-7);
781	f.SetRange(150, 1500);
782	f.SetLineColor(kBlack);
783	f.SetLineStyle(kDashed);
784	f.DrawCopy("same");
785
786	// Plot other spectra from HEGRA
787	f.SetParameter(0, -2.61);
788	f.SetParameter(1, 2.7e-7);
789	f.SetRange(1000, 20000);
790	f.SetLineColor(kBlack);
791	f.SetLineStyle(kDashed);
792	f.DrawCopy("same");
793	*/
794	}
795
796	// --------------------------------------------------------------------------
797	//
798	// Scale some image parameter plots using the scale factor and plot them
799	// together with the corresponding MC histograms.
800	// Called from DisplaySize
801	//
802	Bool_t MJSpectrum::PlotSame(MStatusArray &arr, MParList &plist, const char name, const char tab, const char *plot, Double_t scale) const
803	{
804	TString same(name);
805	same += "Same";
806
807	TH1 h1 = (TH1)arr.FindObjectInCanvas(name, "TH1F", tab);
808	TH1 h2 = (TH1)arr.FindObjectInCanvas(same, "TH1F", tab);
809	if (!h1 \|\| !h2)
810	return kFALSE;
811
812	TObject *obj = plist.FindObject(plot);
813	if (!obj)
814	{
815	*fLog << warn << plot << " not in parameter list... skipping." << endl;
816	return kFALSE;
817	}
818
819	TH1 h3 = (TH1)obj->FindObject(name);
820	if (!h3)
821	{
822	*fLog << warn << name << " not found in " << plot << "... skipping." << endl;
823	return kFALSE;
824	}
825
826
827	const MAlphaFitter fit = (MAlphaFitter)plist.FindObject("MAlphaFitter");
828	const Double_t ascale = fit ? fit->GetScaleFactor() : 1;
829
830	gPad->SetBorderMode(0);
831	h2->SetLineColor(kBlack);
832	h3->SetLineColor(kBlue);
833	h2->Add(h1, -ascale);
834
835	//h2->Scale(1./ontime); //h2->Integral());
836	h3->Scale(scale); //h3->Integral());
837
838	h2->SetMaximum(1.05*TMath::Max(h2->GetMaximum(), h3->GetMaximum()));
839
840	h2 = h2->DrawCopy();
841	h3 = h3->DrawCopy("same");
842
843	// Don't do this on the original object!
844	h2->SetStats(kFALSE);
845	h3->SetStats(kFALSE);
846
847	return kTRUE;
848	}
849
850	// --------------------------------------------------------------------------
851	//
852	// Take a lot of histograms and plot them together in one plot.
853	// Calls PlotSame
854	//
855	Bool_t MJSpectrum::DisplaySize(MParList &plist, Double_t scale) const
856	{
857	*fLog << inf << "Reading from file: " << fPathIn << endl;
858
859	TFile file(fPathIn, "READ");
860	if (!file.IsOpen())
861	{
862	*fLog << err << dbginf << "ERROR - Could not open file " << fPathIn << endl;
863	return kFALSE;
864	}
865
866	file.cd();
867	MStatusArray arr;
868	if (arr.Read()<=0)
869	{
870	*fLog << "MStatusDisplay not found in file... abort." << endl;
871	return kFALSE;
872	}
873
874	TH1 excess = (TH1D)arr.FindObjectInCanvas("Excess", "TH1D", "Hist");
875	if (!excess)
876	return kFALSE;
877
878	// ------------------- Plot excess versus size -------------------
879
880	TCanvas &c = fDisplay->AddTab("Excess");
881	c.Divide(3,2);
882	c.cd(1);
883	gPad->SetBorderMode(0);
884	gPad->SetLogx();
885	gPad->SetLogy();
886	gPad->SetGridx();
887	gPad->SetGridy();
888
889	excess->SetTitle("Excess events vs Size (data/black, mc/blue)");
890	excess = excess->DrawCopy("E2");
891	// Don't do this on the original object!
892	excess->SetStats(kFALSE);
893	excess->SetMarkerStyle(kFullDotMedium);
894	excess->SetFillColor(kBlack);
895	excess->SetFillStyle(0);
896	excess->SetName("Excess ");
897	excess->SetDirectory(0);
898
899	TObject *o=0;
900	if ((o=plist.FindObject("ExcessMC")))
901	{
902	TH1 histsel = (TH1F)o->FindObject("");
903	if (histsel)
904	{
905	if (scale<0)
906	scale = excess->Integral()/histsel->Integral();
907
908	histsel->Scale(scale);
909	histsel->SetLineColor(kBlue);
910	histsel->SetBit(kCanDelete);
911	histsel = histsel->DrawCopy("E1 same");
912	// Don't do this on the original object!
913	histsel->SetStats(kFALSE);
914
915	fLog->Separator("Kolmogorov Test");
916	histsel->KolmogorovTest(excess, "DX");
917	fLog->Separator("Chi^2 Test");
918	const Double_t p = histsel->Chi2Test(excess, "P");
919
920	TLatex tex;
921	tex.SetBit(TLatex::kTextNDC);
922	tex.DrawLatex(0.75, 0.93, Form("P(\\chi^{2})=%.0f%%", p*100));
923	}
924	}
925
926	// -------------- Comparison of Image Parameters --------------
927	c.cd(2);
928	PlotSame(arr, plist, "Dist", "HilSrc", "MHHilSrcMCPost", scale);
929
930	c.cd(3);
931	PlotSame(arr, plist, "Length", "PostCut", "MHHillasMCPost", scale);
932
933	c.cd(4);
934	PlotSame(arr, plist, "M3l", "HilExt", "MHHilExtMCPost", scale);
935
936	c.cd(5);
937	PlotSame(arr, plist, "Conc1", "NewPar", "MHNewParMCPost", scale);
938
939	c.cd(6);
940	PlotSame(arr, plist, "Width", "PostCut", "MHHillasMCPost", scale);
941
942	return kTRUE;
943	}
944
945	Bool_t MJSpectrum::Process(const MDataSet &set)
946	{
947	if (!set.IsValid())
948	{
949	*fLog << err << "ERROR - DataSet invalid!" << endl;
950	return kFALSE;
951	}
952
953	CheckEnv();
954
955	// --------------------------------------------------------------------------------
956
957	*fLog << inf;
958	fLog->Separator(GetDescriptor());
959	*fLog << "Compile Monte Carlo Sample (data set " << set.GetName() << ")" << endl;
960	*fLog << endl;
961
962	MBinning bins1("BinningAlpha");
963	MBinning bins2("BinningEnergyEst");
964	MBinning bins3("BinningTheta");
965	MBinning bins4("BinningFalseSource");
966	MBinning bins5("BinningWidth");
967	MBinning bins6("BinningLength");
968	MBinning bins7("BinningDist");
969	MBinning bins8("BinningMaxDist");
970	MBinning bins9("BinningM3Long");
971	MBinning bins0("BinningConc1");
972
973	MAlphaFitter fit;
974
975	MParList plist;
976	plist.AddToList(&bins1);
977	plist.AddToList(&bins3);
978	plist.AddToList(&bins4);
979	plist.AddToList(&bins5);
980	plist.AddToList(&bins6);
981	plist.AddToList(&bins7);
982	plist.AddToList(&bins8);
983	plist.AddToList(&bins9);
984	plist.AddToList(&bins0);
985	plist.AddToList(&fit);
986
987	TH1D temp1, size;
988	const Float_t ontime = ReadInput(plist, temp1, size);
989	if (ontime<0)
990	{
991	*fLog << err << GetDescriptor() << ": Could not determin effective on time..." << endl;
992	return kFALSE;
993	}
994
995	plist.AddToList(&bins2);
996
997	MMcSpectrumWeight weight;
998	if (!InitWeighting(set, weight))
999	return kFALSE;
1000
1001	PrintSetup(fit);
1002	bins3.SetEdges(temp1, 'x');
1003
1004	TH1D temp2(temp1);
1005	if (!ReadOrigMCDistribution(set, temp2, weight))
1006	return kFALSE;
1007
1008	if (!GetThetaDistribution(temp1, temp2))
1009	return kFALSE;
1010
1011	if (!fNoThetaWeights)
1012	weight.SetWeightsZd(&temp1);
1013
1014	TH1D excess;
1015	if (!Refill(plist, excess))
1016	return kFALSE;
1017
1018	TH2D hist;
1019	MH3 mh1("MMcEvtBasic.fTelescopeTheta*kRad2Deg", "MMcEvtBasic.fEnergy");
1020	if (fSimpleMode)
1021	{
1022	hist.UseCurrentStyle();
1023	MH::SetBinning(&hist, &bins3/temp1.GetXaxis()/, &bins2/excess.GetXaxis()/);
1024	if (!ReadOrigMCDistribution(set, hist, weight))
1025	return kFALSE;
1026
1027	if (!fRawMc)
1028	{
1029	for (int y=0; y<hist.GetNbinsY(); y++)
1030	for (int x=0; x<hist.GetNbinsX(); x++)
1031	hist.SetBinContent(x, y, hist.GetBinContent(x, y)*temp1.GetBinContent(x));
1032	//hist.SetEntries(hist.Integral());
1033	}
1034	}
1035	else
1036	{
1037	weight.SetNameMcEvt("MMcEvtBasic");
1038	if (!IntermediateLoop(plist, mh1, temp1, set, weight))
1039	return kFALSE;
1040	weight.SetNameMcEvt();
1041	}
1042
1043	DisplayResult(fSimpleMode ? hist : (TH2D&)mh1.GetHist());
1044
1045	// ------------------------- Final loop --------------------------
1046
1047	*fLog << endl;
1048	fLog->Separator("Calculate efficiencies");
1049	*fLog << endl;
1050
1051	MTaskList tlist2;
1052	plist.Replace(&tlist2);
1053
1054	MReadMarsFile read("Events");
1055	read.DisableAutoScheme();
1056	set.AddFilesOn(read);
1057
1058	// Selector to get correct (final) theta-distribution
1059	temp1.SetXTitle("MPointingPos.fZd");
1060
1061	MH3 mh3(temp1);
1062
1063	MFEventSelector2 sel2(mh3);
1064	sel2.SetHistIsProbability();
1065
1066	MContinue contsel(&sel2);
1067	contsel.SetInverted();
1068
1069	// Get correct source position
1070	//MSrcPosCalc calc;
1071
1072	// Calculate corresponding Hillas parameters
1073	/*
1074	MHillasCalc hcalc1;
1075	MHillasCalc hcalc2("MHillasCalcAnti");
1076	hcalc1.SetFlags(MHillasCalc::kCalcHillasSrc);
1077	hcalc2.SetFlags(MHillasCalc::kCalcHillasSrc);
1078	hcalc2.SetNameHillasSrc("MHillasSrcAnti");
1079	hcalc2.SetNameSrcPosCam("MSrcPosAnti");
1080	*/
1081
1082	// Fill collection area and energy estimator (unfolding)
1083	// Make sure to use the same binning for MHCollectionArea and MHEnergyEst
1084	MHCollectionArea area;
1085	area.SetHistAll(fSimpleMode ? hist : (TH2D&)mh1.GetHist());
1086	MHEnergyEst hest;
1087
1088	MFillH fill3(&area, "", "FillCollectionArea");
1089	MFillH fill4(&hest, "", "FillEnergyEst");
1090	MFillH fill5("MHThreshold", "", "FillThreshold");
1091	fill3.SetWeight();
1092	fill4.SetWeight();
1093	fill5.SetWeight();
1094	fill3.SetNameTab("ColArea");
1095	fill4.SetNameTab("E-Est");
1096	fill5.SetNameTab("Threshold");
1097
1098	MH3 hsize("MHillas.fSize");
1099	hsize.SetName("ExcessMC");
1100	hsize.Sumw2();
1101
1102	MBinning bins(size, "BinningExcessMC");
1103	plist.AddToList(&hsize);
1104	plist.AddToList(&bins);
1105
1106	MFillH fill1a("MHHillasMCPre [MHHillas]", "MHillas", "FillHillasPre");
1107	MFillH fill2a("MHHillasMCPost [MHHillas]", "MHillas", "FillHillasPost");
1108	MFillH fill3a("MHVsSizeMCPost [MHVsSize]", "MHillasSrc", "FillVsSizePost");
1109	MFillH fill4a("MHHilExtMCPost [MHHillasExt]", "MHillasSrc", "FillHilExtPost");
1110	MFillH fill5a("MHHilSrcMCPost [MHHillasSrc]", "MHillasSrc", "FillHilSrcPost");
1111	MFillH fill6a("MHImgParMCPost [MHImagePar]", "MImagePar", "FillImgParPost");
1112	MFillH fill7a("MHNewParMCPost [MHNewImagePar]", "MNewImagePar", "FillNewParPost");
1113	MFillH fill8a("ExcessMC [MH3]", "", "FillExcessMC");
1114	fill1a.SetNameTab("PreCut");
1115	fill2a.SetNameTab("PostCut");
1116	fill3a.SetNameTab("VsSize");
1117	fill4a.SetNameTab("HilExt");
1118	fill5a.SetNameTab("HilSrc");
1119	fill6a.SetNameTab("ImgPar");
1120	fill7a.SetNameTab("NewPar");
1121	fill8a.SetBit(MFillH::kDoNotDisplay);
1122	fill1a.SetWeight();
1123	fill2a.SetWeight();
1124	fill3a.SetWeight();
1125	fill4a.SetWeight();
1126	fill5a.SetWeight();
1127	fill6a.SetWeight();
1128	fill7a.SetWeight();
1129	fill8a.SetWeight();
1130
1131	MEnergyEstimate est;
1132	MTaskEnv taskenv1("EstimateEnergy");
1133	taskenv1.SetDefault(fEstimateEnergy ? fEstimateEnergy : &est);
1134
1135	tlist2.AddToList(&read);
1136	if (!fRawMc && fNoThetaWeights)
1137	tlist2.AddToList(&contsel);
1138	//tlist2.AddToList(&calc);
1139	//tlist2.AddToList(&hcalc1);
1140	//tlist2.AddToList(&hcalc2);
1141	tlist2.AddToList(&weight);
1142	tlist2.AddToList(&fill1a);
1143	tlist2.AddToList(fCut0);
1144	tlist2.AddToList(fCut1);
1145	tlist2.AddToList(fCut2);
1146	tlist2.AddToList(fCut3);
1147	tlist2.AddToList(&taskenv1);
1148	tlist2.AddToList(&fill3);
1149	tlist2.AddToList(&fill4);
1150	tlist2.AddToList(&fill5);
1151	tlist2.AddToList(&fill2a);
1152	tlist2.AddToList(&fill3a);
1153	tlist2.AddToList(&fill4a);
1154	tlist2.AddToList(&fill5a);
1155	tlist2.AddToList(&fill6a);
1156	tlist2.AddToList(&fill7a);
1157	tlist2.AddToList(&fill8a);
1158	//tlist2.AddToList(&fill9a);
1159
1160	MEvtLoop loop2("FillMonteCarlo"); // *** fName ***
1161	loop2.SetParList(&plist);
1162	loop2.SetDisplay(fDisplay);
1163	loop2.SetLogStream(fLog);
1164
1165	if (!SetupEnv(loop2))
1166	return kFALSE;
1167
1168	if (!loop2.Eventloop(fMaxEvents))
1169	{
1170	*fLog << err << GetDescriptor() << ": Processing of MC-data failed." << endl;
1171	return kFALSE;
1172	}
1173
1174	if (!loop2.GetDisplay())
1175	{
1176	*fLog << err << GetDescriptor() << ": Execution stopped by user." << endl;
1177	return kFALSE;
1178	}
1179
1180	gLog.Separator("Energy Estimator");
1181	if (plist.FindObject("EstimateEnergy"))
1182	plist.FindObject("EstimateEnergy")->Print();
1183
1184	gLog.Separator("Spectrum");
1185
1186	// -------------------------- Spectrum ----------------------------
1187
1188	// Calculate and display spectrum (N/TeVsm^2 at 1TeV)
1189	TArrayD res(DisplaySpectrum(area, excess, hest, ontime));
1190
1191	// Spectrum fitted (convert res[1] from TeV to GeV)
1192	TF1 flx("flx", Form("%e*pow(x/1000, %f)", res[1]/1000, res[0]));
1193
1194	// Number of events this spectrum would produce per s and m^2
1195	Double_t n = flx.Integral(weight.GetEnergyMin(), weight.GetEnergyMax());
1196
1197	// scale with effective collection area to get the event rate (N/s)
1198	// scale with the effective observation time to absolute observed events
1199	n = area.GetCollectionAreaAbs()ontime; // N
1200
1201	// Now calculate the scale factor from the number of events
1202	// produced and the number of events which should have been
1203	// observed with our telescope in the time ontime
1204	const Double_t scale = n/area.GetEntries();
1205
1206	// Print normalization constant
1207	cout << "MC normalization factor: " << scale << endl;
1208
1209	// Overlay normalized plots
1210	DisplaySize(plist, scale);
1211
1212	// check if output should be written
1213	if (!fPathOut.IsNull())
1214	fDisplay->SaveAsRoot(fPathOut);
1215
1216	return kTRUE;
1217	}

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

Download in other formats: