Context Navigation

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

Visit:

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

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

Download in other formats: