Context Navigation

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

Visit:

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

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

Download in other formats: