Context Navigation

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

Visit:

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

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

Download in other formats: