Context Navigation

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

Visit:

Last change on this file since 8047 was 8047, checked in by tbretz, 18 years ago
* empty log message *
File size: 35.6 KB

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

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

Download in other formats: