Context Navigation

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

Visit:

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

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

Download in other formats: