Context Navigation

source: releases/Mars.2014.05.26/mjobs/MJSpectrum.cc@ 18459

Visit:

Last change on this file since 18459 was 9851, checked in by tbretz, 14 years ago
Changed MH::SetBinning and similar functions to take references instead of pointers as arguments. For convenience wrappers for the old style functions are provided.
File size: 60.3 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-2008
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 <TGraph.h>
41	#include <TChain.h>
42	#include <TLatex.h>
43	#include <TCanvas.h>
44	#include <TObjArray.h>
45
46	// Environment
47	#include "MLog.h"
48	#include "MLogManip.h"
49
50	#include "MMath.h"
51	#include "MString.h"
52	#include "MDirIter.h"
53
54	#include "MStatusArray.h"
55	#include "MStatusDisplay.h"
56
57	// Container
58	#include "MH3.h"
59	#include "MHn.h"
60	#include "MEnv.h"
61	#include "MBinning.h"
62	#include "MParameters.h"
63	#include "MDataSet.h"
64	#include "MMcCorsikaRunHeader.h"
65
66	// Spectrum
67	#include "MAlphaFitter.h"
68	#include "MHAlpha.h"
69	#include "MHCollectionArea.h"
70	#include "MHEnergyEst.h"
71	#include "MMcSpectrumWeight.h"
72	#include "MHEffectiveOnTime.h"
73
74	// Eventloop
75	#include "MEvtLoop.h"
76	#include "MTaskList.h"
77	#include "MParList.h"
78
79	// Tasks/Filter
80	#include "MReadMarsFile.h"
81	#include "MEnergyEstimate.h"
82	#include "MTaskEnv.h"
83	#include "MFillH.h"
84	#include "MFDataPhrase.h"
85	#include "MFDataMember.h"
86	#include "MContinue.h"
87	#include "MWriteRootFile.h"
88
89	ClassImp(MJSpectrum);
90
91	using namespace std;
92
93	MJSpectrum::MJSpectrum(const char name, const char title)
94	: fCutQ(0), fCut0(0),fCut1(0), fCut2(0), fCut3(0), fCutS(0),
95	fEstimateEnergy(0), fCalcHadronness(0), fCalcDisp(0),
96	fForceTheta(kFALSE), fForceRunTime(kFALSE), fForceOnTimeFit(kFALSE)
97	{
98	fName = name ? name : "MJSpectrum";
99	fTitle = title ? title : "Standard program to calculate spectrum";
100
101	// Make sure that fDisplay is maintained properly
102	// (i.e. removed via RecursiveRemove if closed)
103	gROOT->GetListOfCleanups()->Add(this);
104	SetBit(kMustCleanup);
105	}
106
107	// --------------------------------------------------------------------------
108	//
109	// Delete all the fCut* data members and fCalcHadronness/fEstimateEnergy
110	//
111	MJSpectrum::~MJSpectrum()
112	{
113	if (fCutQ)
114	delete fCutQ;
115	if (fCut0)
116	delete fCut0;
117	if (fCut1)
118	delete fCut1;
119	if (fCut2)
120	delete fCut2;
121	if (fCut3)
122	delete fCut3;
123	if (fCutS)
124	delete fCutS;
125	if (fEstimateEnergy)
126	delete fEstimateEnergy;
127	if (fCalcHadronness)
128	delete fCalcHadronness;
129	if (fCalcDisp)
130	delete fCalcDisp;
131	}
132
133	// --------------------------------------------------------------------------
134	//
135	// Setup a task estimating the energy. The given task is cloned.
136	//
137	void MJSpectrum::SetEnergyEstimator(const MTask *task)
138	{
139	if (fEstimateEnergy)
140	delete fEstimateEnergy;
141	fEstimateEnergy = task ? (MTask*)task->Clone() : 0;
142	}
143
144	// --------------------------------------------------------------------------
145	//
146	// Read a MTask named name into task from the open file. If this task is
147	// required mustexist can be set. Depending on success kTRUE or kFALSE is
148	// returned. If the task is a MContinue SetAllowEmpty is called.
149	// The name of the task is set to name.
150	//
151	Bool_t MJSpectrum::ReadTask(MTask* &task, const char *name, Bool_t mustexist) const
152	{
153	// If a task is set delete task
154	if (task)
155	{
156	delete task;
157	task = 0;
158	}
159
160	// Try to get task from file
161	task = (MTask*)gFile->Get(name);
162	if (!task)
163	{
164	if (!mustexist)
165	return kTRUE;
166	*fLog << err << dbginf << "ERROR - " << name << " doen't exist in file!" << endl;
167	return kFALSE;
168	}
169
170	// Check if task inherits from MTask
171	if (!task->InheritsFrom(MTask::Class()))
172	{
173	*fLog << err << dbginf << "ERROR - " << name << " read doesn't inherit from MTask!" << endl;
174	delete task;
175	return kFALSE;
176	}
177
178	// Set name of task
179	task->SetName(name);
180
181	// SetAllowEmpty for MContinue tasks
182	if (dynamic_cast<MContinue*>(task))
183	dynamic_cast<MContinue*>(task)->SetAllowEmpty();
184
185	return kTRUE;
186	}
187
188	// --------------------------------------------------------------------------
189	//
190	// Print setup used for the MC processing, namely MAlphaFitter ans all fCut*.
191	//
192	void MJSpectrum::PrintSetup(const MAlphaFitter &fit) const
193	{
194	fLog->Separator("Alpha Fitter");
195	*fLog << all;
196	fit.Print("result");
197
198	fLog->Separator("Used Cuts");
199	fCutQ->Print();
200	fCut0->Print();
201	fCut1->Print();
202	fCutS->Print();
203	fCut2->Print();
204	fCut3->Print();
205	}
206
207	// --------------------------------------------------------------------------
208	//
209	// Read the first MMcCorsikaRunHeader from the RunHeaders tree in
210	// the dataset.
211	// The simulated energy range and spectral slope is initialized from
212	// there.
213	// In the following eventloops the forced check in MMcSpectrumWeight
214	// ensures, that the spectral slope and energy range doesn't change.
215	//
216	Bool_t MJSpectrum::InitWeighting(const MDataSet &set, MMcSpectrumWeight &w) const
217	{
218	fLog->Separator("Initialize energy weighting");
219
220	// Read Resources
221	if (!CheckEnv(w))
222	{
223	*fLog << err << "ERROR - Reading resources for MMcSpectrumWeight failed." << endl;
224	return kFALSE;
225	}
226
227	w.Print("new");
228
229	return kTRUE;
230	}
231
232	Float_t MJSpectrum::ReadInput(MParList &plist, TH1D &h1, TH1D &h2)
233	{
234	*fLog << inf << "Reading from file: " << fPathIn << endl;
235
236	TFile file(fPathIn, "READ");
237	if (!file.IsOpen())
238	{
239	*fLog << err << dbginf << "ERROR - Could not open file " << fPathIn << endl;
240	return -1;
241	}
242
243	MStatusArray arr;
244	if (arr.Read()<=0)
245	{
246	*fLog << err << dbginf << "ERROR - MStatusDisplay not found in file... abort." << endl;
247	return -1;
248	}
249
250	//arr.EnableTH1Workaround();
251
252	TH1D vstime = (TH1D)arr.FindObjectInCanvas("Theta", "TH1D", "OnTime");
253	TH1D size = (TH1D)arr.FindObjectInCanvas("Excess", "TH1D", "Hist");
254	if (!vstime)
255	{
256	*fLog << err << dbginf << "ERROR - Theta [TH1D] not found in OnTime-tab... abort." << endl;
257	return -1;
258	}
259	if (!size)
260	{
261	*fLog << err << dbginf << "ERROR - Excess [TH1D] not found in Hist-tab... abort." << endl;
262	return -1;
263	}
264
265	vstime->Copy(h1);
266	size->Copy(h2);
267	h1.SetDirectory(0);
268	h2.SetDirectory(0);
269
270	if (fDisplay)
271	arr.DisplayIn(*fDisplay, "Hist");
272
273	if (!ReadTask(fCutQ, "CutQ"))
274	return -1;
275	if (!ReadTask(fCut0, "Cut0"))
276	return -1;
277	if (!ReadTask(fCut1, "Cut1"))
278	return -1;
279	if (!ReadTask(fCut2, "Cut2"))
280	return -1;
281	if (!ReadTask(fCut3, "Cut3"))
282	return -1;
283	if (!ReadTask(fCalcHadronness, "CalcHadronness", kFALSE))
284	return -1;
285	if (!ReadTask(fCalcDisp, "CalcDisp", kFALSE))
286	return -1;
287
288	TObjArray arrread;
289
290	TIter Next(plist);
291	TObject *o=0;
292	while ((o=Next()))
293	if (o->InheritsFrom(MBinning::Class()))
294	arrread.Add(o);
295
296	arrread.Add(plist.FindObject("MAlphaFitter"));
297
298	if (!ReadContainer(arrread))
299	return -1;
300
301	// We don't have to proceed. We later overwrite the result anyway
302	if (fForceOnTimeFit)
303	return 0;
304
305	const Double_t ufl = vstime->GetBinContent(0);
306	const Double_t ofl = vstime->GetBinContent(vstime->GetNbinsX()+1);
307	const Double_t eff = vstime->Integral()+ufl+ofl;
308	if (ufl>0)
309	{
310	if (vstime->GetBinLowEdge(1)<=0)
311	{
312	*fLog << err << "ERROR - Undeflow bin of OnTime histogram not empty as it ought to be." << endl;
313	return -1;
314	}
315	fLog << warn << "WARNING - " << Form("%.1f%% (%.0fs)", 100ufl/eff, ufl) << " of the eff. observation time is in the underflow bin." << endl;
316	}
317	if (ofl>0)
318	fLog << warn << "WARNING - " << Form("%.1f%% (%.0fs)", 100ofl/eff, ofl) << " of the eff. observation time is in the overflow bin." << endl;
319
320	if (!fForceRunTime)
321	return eff;
322
323	TH1D time = (TH1D)arr.FindObjectInCanvas("ExcessTime", "TH1D", "Hist");
324	if (!time)
325	{
326	*fLog << err;
327	*fLog << "ERROR - ExcessTime [TH1D] not found in Hist-tab... abort." << endl;
328	*fLog << " Did you try to process Monte Carlos with --force-runtime?" <<endl;
329	return -1;
330	}
331
332	const Double_t obs = time->GetXaxis()->GetXmax()-time->GetXaxis()->GetXmin();
333
334	*fLog << inf;
335	*fLog << "Total run time: " << obs/60 << "min" << endl;
336	fLog << "Eff. obs. time: " << eff/60 << "min (" << Form("%.1f%%", 100eff/obs) << ")" << endl;
337
338	return obs;
339	}
340
341	// --------------------------------------------------------------------------
342	//
343	// return maximum value of MMcRunHeader.fImpactMax stored in the RunHeaders
344	// of the files from the MC dataset
345	//
346	Bool_t MJSpectrum::AnalyzeMC(const MDataSet &set, Float_t &impactmax, Float_t &emin/, Float_t emax/) const
347	{
348	if (fDisplay)
349	fDisplay->SetStatusLine1("Analyzing Monte Carlo headers...");
350
351	// ----- Create chain -----
352	*fLog << inf << "Getting files... " << flush;
353	TChain chain("RunHeaders");
354	if (!set.AddFilesOn(chain))
355	return kFALSE;
356	*fLog << "done. " << endl;
357
358	*fLog << all;
359	*fLog << "Searching for maximum impact... " << flush;
360	impactmax = chain.GetMaximum("MMcRunHeader.fImpactMax");
361	*fLog << "found " << impactmax/100 << "m" << endl;
362
363	*fLog << "Searching for minimum lower energy limit... " << flush;
364	emin = chain.GetMinimum("MMcCorsikaRunHeader.fELowLim");
365	*fLog << "found " << emin << "GeV" << endl;
366
367	*fLog << "Searching for maximum lower energy limit... " << flush;
368	const Float_t emin2 = chain.GetMaximum("MMcCorsikaRunHeader.fELowLim");
369	*fLog << "found " << emin2 << "GeV" << endl;
370
371	// Need a check for the upper energy LIMIT?!?
372
373	if (emin2>emin)
374	{
375	*fLog << warn;
376	*fLog << "WARNING - You are using files with different lower limits for the simulated" << endl;
377	*fLog << " energy, i.e. that the collection area and your correction" << endl;
378	*fLog << " coefficients between " << emin << "GeV and ";
379	*fLog << emin2 << "GeV might be wrong." << endl;
380	}
381
382	/*
383	*fLog << "Getting maximum energy... " << flush;
384	emax = chain.GetMaximum("MMcCorsikaRunHeader.fEUppLim");
385	*fLog << "found " << emax << "GeV" << endl;
386	*/
387	return kTRUE;
388	}
389
390	Bool_t MJSpectrum::ReadOrigMCDistribution(const MDataSet &set, TH1 &h, MMcSpectrumWeight &weight) const
391	{
392	// Some debug output
393	*fLog << all << endl;
394	fLog->Separator("Compiling original MC distribution");
395
396	Float_t impactmax=0, Emin=0;
397	if (!AnalyzeMC(set, impactmax, Emin))
398	return kFALSE;
399
400	*fLog << inf << "Getting files... " << flush;
401	MDirIter iter;
402	if (!set.AddFilesOn(iter))
403	return kFALSE;
404	*fLog << "done. " << endl;
405
406	const Int_t tot = iter.GetNumEntries();
407
408	// Prepare histogram
409	h.Reset();
410	h.Sumw2();
411	if (h.InheritsFrom(TH2::Class()))
412	{
413	h.SetNameTitle("ThetaEMC", "Event-Distribution vs Theta and Energy for MC (produced)");
414	h.SetXTitle("\\Theta [\\circ]");
415	h.SetYTitle("E [GeV]");
416	h.SetZTitle(MString::Format("Counts normalized to r=%.0fm", impactmax/100));
417	}
418	else
419	{
420	h.SetNameTitle("ThetaMC", "Event-Distribution vs Theta for MC (produced)");
421	h.SetXTitle("\\Theta [\\circ]");
422	h.SetYTitle(MString::Format("Counts normalized to r=%.0fm", impactmax/100));
423	}
424	h.SetDirectory(0);
425
426	const TString cont = h.InheritsFrom(TH2::Class()) ?
427	"MMcEvtBasic.fEnergy:MMcEvtBasic.fTelescopeTheta*TMath::RadToDeg()>>ThetaEMC" :
428	"MMcEvtBasic.fTelescopeTheta*TMath::RadToDeg()>>ThetaMC";
429
430	if (fDisplay)
431	fDisplay->SetStatusLine1("Reading OriginalMC distribution...");
432
433	TH1 hfill = (TH1)h.Clone(h.InheritsFrom(TH2::Class())?"ThetaEMC":"ThetaMC");
434	hfill->SetDirectory(0);
435
436	*fLog << all << "Compiling simulated source spectrum... please stand by, this may take a while." << endl;
437
438	Double_t evts = 0;
439	Int_t num = 0;
440
441	// Reading this with a eventloop is five times slower :(
442	TString fname;
443	while (fDisplay)
444	{
445	if (fDisplay)
446	fDisplay->SetProgressBarPosition(Float_t(num++)/tot);
447
448	// Get next filename
449	fname = iter.Next();
450	if (fname.IsNull())
451	break;
452
453	if (fDisplay)
454	fDisplay->SetStatusLine2(fname);
455
456	// open file
457	TFile file(fname);
458	if (file.IsZombie())
459	{
460	*fLog << err << "ERROR - Couldn't open file " << fname << endl;
461	return kFALSE;
462	}
463
464	// Get tree OriginalMC
465	TTree t = dynamic_cast<TTree>(file.Get("OriginalMC"));
466	if (!t)
467	{
468	*fLog << err << "ERROR - File " << fname << " doesn't contain tree OriginalMC." << endl;
469	return kFALSE;
470	}
471	if (t->GetEntries()==0)
472	{
473	*fLog << warn << "WARNING - OriginalMC of " << fname << " empty." << endl;
474	continue;
475	}
476
477	// Get tree RunHeaders
478	TTree rh = dynamic_cast<TTree>(file.Get("RunHeaders"));
479	if (!rh)
480	{
481	*fLog << err << "ERROR - File " << fname << " doesn't contain tree RunHeaders." << endl;
482	return kFALSE;
483	}
484	if (rh->GetEntries()!=1)
485	{
486	*fLog << err << "ERROR - RunHeaders of " << fname << " doesn't contain exactly one entry." << endl;
487	return kFALSE;
488	}
489
490	// Get corsika run header
491	MMcCorsikaRunHeader *head=0;
492	rh->SetBranchAddress("MMcCorsikaRunHeader.", &head);
493	rh->GetEntry(0);
494	if (!head)
495	{
496	*fLog << err << "ERROR - Couldn't read MMcCorsikaRunHeader from " << fname << "." << endl;
497	return kFALSE;
498	}
499
500	// Get the maximum impact parameter of this file. Due to different
501	// production areas an additional scale-factor is applied.
502	//
503	// Because it is assumed that the efficiency outside the production
504	// area is nearly zero no additional weight has to be applied to the
505	// events after cuts. For the events before cuts it is fair to use
506	// weights... maybe filling the residual impact with unweighted
507	// events would be better?!? (Not that the weighting might be
508	// less correct with low statistics, because it could pronounce
509	// a fluctuation)
510	const Double_t impact = rh->GetMaximum("MMcRunHeader.fImpactMax");
511	const Double_t scale = impactmax/impact;
512
513	// Propagate the run header to MMcSpectrumWeight
514	if (!weight.Set(*head))
515	{
516	*fLog << err << "ERROR - Initializing MMcSpectrumWeight from " << fname << " failed." << endl;
517	return kFALSE;
518	}
519
520	// Get the corresponding rule for the weights
521	const TString weights(weight.GetFormulaWeights("MMcEvtBasic"));
522
523	// No we found everything... go on reading contents
524	*fLog << inf2 << "Reading OriginalMC of " << fname << endl;
525
526	// Fill histogram from tree
527	hfill->SetDirectory(&file);
528	if (t->Draw(cont, weights, "goff")<0)
529	{
530	*fLog << err << "ERROR - Reading OriginalMC failed." << endl;
531	return kFALSE;
532	}
533	hfill->SetDirectory(0);
534
535	evts += hfill->GetEntries();
536
537	// ----- Complete incomplete energy ranges -----
538	// ----- and apply production area weights -----
539	weight.CompleteEnergySpectrum(hfill, Emin, scalescale);
540
541	// Add weighted data from file to final histogram
542	h.Add(hfill);
543	}
544	delete hfill;
545
546	*fLog << all << "Total number of events from files in Monte Carlo dataset: " << evts << endl;
547	if (fDisplay)
548	fDisplay->SetStatusLine2("done.");
549
550	if (!fDisplay \|\| h.GetEntries()>0)
551	return kTRUE;
552
553	*fLog << err << "ERROR - Histogram with distribution from OriginalMC empty..." << endl;
554	return kFALSE;
555	}
556
557	void MJSpectrum::GetThetaDistribution(TH1D &temp1, TH2D &hist2) const
558	{
559	TH1D &temp2 = *hist2.ProjectionX();
560
561	// Display some stuff
562	if (fDisplay)
563	{
564	TCanvas &c = fDisplay->AddTab("ZdDist");
565	c.Divide(2,2);
566
567	// On-Time vs. Theta
568	c.cd(1);
569	gPad->SetBorderMode(0);
570	gPad->SetGridx();
571	gPad->SetGridy();
572	temp1.DrawCopy();
573
574	// Number of MC events (produced) vs Theta
575	c.cd(2);
576	gPad->SetBorderMode(0);
577	gPad->SetGridx();
578	gPad->SetGridy();
579	temp2.SetName("NVsTheta");
580	temp2.DrawCopy("hist");
581
582	c.cd(4);
583	gPad->SetBorderMode(0);
584	gPad->SetGridx();
585	gPad->SetGridy();
586
587	c.cd(3);
588	gPad->SetBorderMode(0);
589	gPad->SetGridx();
590	gPad->SetGridy();
591	}
592
593	// Calculate the Probability
594	temp1.Divide(&temp2);
595	temp1.Scale(1./temp1.Integral(1, temp1.GetNbinsX()+1));
596
597	// Some cosmetics: Name, Axis, etc.
598	temp1.SetName("ProbVsTheta");
599	temp1.SetTitle("Probability vs. Zenith Angle to choose MC events");
600	temp1.SetYTitle("Probability");
601	if (fDisplay)
602	temp1.DrawCopy("hist");
603
604	delete &temp2;
605	}
606
607	// --------------------------------------------------------------------------
608	//
609	// Display the final theta distribution.
610	//
611	Bool_t MJSpectrum::DisplayResult(const TH2D &h2) const
612	{
613	if (!fDisplay \|\| !fDisplay->CdCanvas("ZdDist"))
614	{
615	*fLog << err << "ERROR - Display or tab ZdDist vanished... abort." << endl;
616	return kFALSE;
617	}
618
619	TH1D &proj = *h2.ProjectionX();
620	proj.SetNameTitle("ThetaFinal", "Final Theta Distribution");
621	proj.SetXTitle("\\Theta [\\circ]");
622	proj.SetYTitle("Counts");
623	proj.SetLineColor(kBlue);
624	proj.SetDirectory(0);
625	proj.SetBit(kCanDelete);
626
627	TVirtualPad *pad = gPad;
628
629	pad->cd(4);
630	proj.DrawCopy();
631
632	pad->cd(1);
633	TH1D theta = (TH1D)gPad->FindObject("Theta");
634	if (!theta)
635	{
636	*fLog << err << "ERROR - Theta-Histogram vanished... cannot proceed." << endl;
637	return kFALSE;
638	}
639
640	// Check whether histogram is empty
641	if (proj.GetMaximum()==0)
642	{
643	*fLog << err;
644	*fLog << "ERROR - The Zenith Angle distribution of your Monte Carlos doesn't overlap" << endl;
645	*fLog << " with the Zenith Angle distribution of your observation." << endl;
646	*fLog << " Maybe the energy binning is undefined or wrong (from ";
647	*fLog << h2.GetYaxis()->GetXmin() << "GeV to " << h2.GetYaxis()->GetXmax() << "GeV)" << endl;
648	theta->SetLineColor(kRed);
649	return kFALSE;;
650	}
651
652	// scale histogram and set new maximum for display
653	proj.Scale(theta->GetMaximum()/proj.GetMaximum());
654	theta->SetMaximum(1.05*TMath::Max(theta->GetMaximum(), proj.GetMaximum()));
655
656	// draw project
657	proj.Draw("same");
658
659	// Compare both histograms
660	*fLog << inf << "Comparing theta distributions for data and MCs." << endl;
661
662	const Double_t prob = proj.Chi2Test(theta, "P");
663	if (prob==1)
664	return kTRUE;
665
666	if (prob>0.99)
667	{
668	*fLog << inf;
669	*fLog << "The Zenith Angle distribution of your Monte Carlos fits well" << endl;
670	*fLog << "with the Zenith Angle distribution of your observation." << endl;
671	*fLog << "The probability for identical Theta distributions is " << prob << endl;
672	return kTRUE;
673	}
674
675	if (prob<0.01)
676	{
677	*fLog << err;
678	*fLog << "ERROR - The Zenith Angle distribution of your Monte Carlos does not fit" << endl;
679	*fLog << " with the Zenith Angle distribution of your observation." << endl;
680	*fLog << " The probability for identical Theta distributions is " << prob << endl;
681	if (!fForceTheta)
682	*fLog << " To force processing use --force-theta (with care!)" << endl;
683	theta->SetLineColor(kRed);
684	return fForceTheta;
685	}
686
687	*fLog << warn;
688	*fLog << "WARNING - The Zenith Angle distribution of your Monte Carlos doesn't fits well" << endl;
689	*fLog << " with the Zenith Angle distribution of your observation." << endl;
690	*fLog << " The probability for identical Theta distributions is " << prob << endl;
691	return kTRUE;
692	}
693
694	TString MJSpectrum::GetHAlpha() const
695	{
696	TString cls("MHAlpha");
697	if (!fDisplay)
698	return cls;
699
700	TCanvas *c = fDisplay->GetCanvas("Hist");
701	if (!c)
702	return cls;
703
704	TIter Next(c->GetListOfPrimitives());
705	TObject *obj=0;
706	while ((obj=Next()))
707	if (obj->InheritsFrom(MHAlpha::Class()))
708	break;
709
710	return obj ? TString(obj->ClassName()) : cls;
711	}
712
713	// --------------------------------------------------------------------------
714	//
715	// Fills the excess histogram (vs E-est) from the events stored in the
716	// ganymed result file and therefor estimates the energy.
717	//
718	// The resulting histogram excess-vs-energy ist copied into h2.
719	//
720	Bool_t MJSpectrum::Refill(MParList &plist, TH1D &h2)/const/
721	{
722	// Now fill the histogram
723	*fLog << endl;
724	fLog->Separator("Refill Excess");
725	*fLog << endl;
726
727	MTaskList tlist;
728	plist.AddToList(&tlist);
729
730	MReadTree read("Events");
731	read.DisableAutoScheme();
732	read.AddFile(fPathIn);
733	/*
734	MTaskEnv taskenv0("CalcDisp");
735	taskenv0.SetDefault(fCalcDisp);
736
737	MTaskEnv taskenv1("CalcHadronness");
738	taskenv1.SetDefault(fCalcHadronness);
739	*/
740	MEnergyEstimate est;
741	MTaskEnv taskenv2("EstimateEnergy");
742	taskenv2.SetDefault(fEstimateEnergy ? fEstimateEnergy : &est);
743
744	MContinue *cont = new MContinue("", "CutS");
745	cont->SetAllowEmpty();
746
747	if (fCutS)
748	delete fCutS;
749	fCutS = cont;
750
751	// Try to find the class used to determine the signal!
752	const TString cls = GetHAlpha();
753
754	// FIXME: Create HistE and HistEOff to be able to modify it from
755	// the resource file.
756	MFillH fill1(Form("HistEOff [%s]", cls.Data()), "", "FillHistEOff");
757	MFillH fill2(Form("HistE [%s]", cls.Data()), "", "FillHistE");
758
759	// Fill a new MHEffectiveOnTime. It can either be used to
760	// re-calculate the on-time or just for manual cross check
761	MFillH fillT("MHEffectiveOnTime", "MTime", "FillOnTime");
762	fillT.SetNameTab("OnTime");
763
764	MFDataMember f0("DataType.fVal", '<', 0.5, "FilterOffData");
765	MFDataMember f1("DataType.fVal", '>', 0.5, "FilterOnData");
766
767	fill1.SetFilter(&f0);
768	fill2.SetFilter(&f1);
769	fillT.SetFilter(&f1);
770
771	tlist.AddToList(&read);
772	//tlist.AddToList(&taskenv0); // not necessary, stored in file!
773	//tlist.AddToList(&taskenv1); // not necessary, stored in file!
774	tlist.AddToList(&f1);
775	tlist.AddToList(&fillT);
776	tlist.AddToList(fCutS);
777	tlist.AddToList(&taskenv2);
778	tlist.AddToList(&f0);
779	tlist.AddToList(&fill1);
780	tlist.AddToList(&fill2);
781
782	// by setting it here it is distributed to all consecutive tasks
783	tlist.SetAccelerator(MTask::kAccDontReset\|MTask::kAccDontTime);
784
785	MEvtLoop loop("RefillExcess"); // *** fName ***
786	loop.SetParList(&plist);
787	loop.SetDisplay(fDisplay);
788	loop.SetLogStream(fLog);
789
790	if (!SetupEnv(loop))
791	return kFALSE;
792
793	if (!loop.Eventloop())
794	{
795	*fLog << err << GetDescriptor() << ": Refilling of data failed." << endl;
796	return kFALSE;
797	}
798
799	if (!loop.GetDisplay())
800	{
801	*fLog << err << GetDescriptor() << ": Execution stopped by user." << endl;
802	return kFALSE;
803	}
804
805	const MHAlpha halpha = (MHAlpha )plist.FindObject("HistE");
806	if (!halpha)
807	{
808	*fLog << err << GetDescriptor() << ": HistE [MHAlpha] not found... abort." << endl;
809	return kFALSE;
810	}
811
812	halpha->GetHEnergy().Copy(h2);
813	h2.SetDirectory(0);
814
815	return kTRUE;
816	}
817
818	/*
819	Bool_t MJSpectrum::IntermediateLoop(MParList &plist, MH3 &mh1, TH1D &temp1, const MDataSet &set, MMcSpectrumWeight &weight) const
820	{
821	MTaskList tlist1;
822	plist.Replace(&tlist1);
823
824	MReadMarsFile readmc("OriginalMC");
825	//readmc.DisableAutoScheme();
826	if (!set.AddFilesOn(readmc))
827	return kFALSE;
828
829	readmc.EnableBranch("MMcEvtBasic.fTelescopeTheta");
830	readmc.EnableBranch("MMcEvtBasic.fEnergy");
831
832	mh1.SetLogy();
833	mh1.SetLogz();
834	mh1.SetName("ThetaE");
835
836	MFillH fill0(&mh1);
837	//fill0.SetDrawOption("projx only");
838
839	MBinning bins2 = (MBinning)plist.FindObject("BinningEnergyEst");
840	MBinning bins3 = (MBinning)plist.FindObject("BinningTheta");
841	if (bins2 && bins3)
842	{
843	bins2->SetName("BinningThetaEY");
844	bins3->SetName("BinningThetaEX");
845	}
846	tlist1.AddToList(&readmc);
847	tlist1.AddToList(&weight);
848
849	temp1.SetXTitle("MMcEvtBasic.fTelescopeTheta*kRad2Deg");
850	MH3 mh3mc(temp1);
851
852	MFEventSelector2 sel1(mh3mc);
853	sel1.SetHistIsProbability();
854
855	fill0.SetFilter(&sel1);
856
857	//if (!fRawMc)
858	tlist1.AddToList(&sel1);
859	tlist1.AddToList(&fill0);
860
861	// by setting it here it is distributed to all consecutive tasks
862	tlist1.SetAccelerator(MTask::kAccDontReset\|MTask::kAccDontTime);
863
864	MEvtLoop loop1("IntermediateLoop"); // *** fName ***
865	loop1.SetParList(&plist);
866	loop1.SetLogStream(fLog);
867	loop1.SetDisplay(fDisplay);
868
869	if (!SetupEnv(loop1))
870	return kFALSE;
871
872	if (!loop1.Eventloop(fMaxEvents))
873	{
874	*fLog << err << GetDescriptor() << ": Processing of MC-data failed." << endl;
875	return kFALSE;
876	}
877
878	if (!loop1.GetDisplay())
879	{
880	*fLog << err << GetDescriptor() << ": Execution stopped by user." << endl;
881	return kFALSE;
882	}
883
884	if (bins2 && bins3)
885	{
886	bins2->SetName("BinningEnergyEst");
887	bins3->SetName("BinningTheta");
888	}
889
890	return kTRUE;
891	}
892	*/
893
894	TString MJSpectrum::FormFloat(Double_t d)
895	{
896	TString s;
897	s += d;
898	return s.Strip(TString::kLeading);
899	}
900
901	TString MJSpectrum::FormFlux(const TF1 &f, const char *unit)
902	{
903	Double_t p0 = -f.GetParameter(0);
904	Double_t p1 = f.GetParameter(1);
905
906	Double_t e0 = f.GetParError(0);
907	Double_t e1 = f.GetParError(1);
908
909	MMath::Format(p0, e0);
910	MMath::Format(p1, e1);
911
912	const Int_t i = TMath::FloorNint(TMath::Log10(p1));
913	const Double_t exp = TMath::Power(10, i);
914
915	TString str = MString::Format("(%s #pm %s)·10^{%d} ",
916	FormFloat(p1/exp).Data(), FormFloat(e1/exp).Data(), i);
917
918	str += MString::Format("#left(#frac{E}{%s}#right)^{-%s #pm %s}", unit,
919	FormFloat(p0).Data(), FormFloat(e0).Data());
920
921	str += " TeV^{-1} m^{-2} s^{-1}";
922
923	return str;
924	}
925
926	TString MJSpectrum::FormString(const TF1 &f, Byte_t type)
927	{
928	switch (type)
929	{
930	case 0:
931	return FormFlux(f, "500GeV");
932	case 1:
933	return MString::Format("\\chi^{2}/NDF=%.2f/%d", f.GetChisquare(),f.GetNDF());
934	case 2:
935	return MString::Format("P=%.0f%%", 100*TMath::Prob(f.GetChisquare(), f.GetNDF()));
936	}
937	return "";
938	}
939
940	TArrayD MJSpectrum::FitSpectrum(TH1D &spectrum) const
941	{
942	Axis_t lo, hi;
943	MH::GetRangeUser(spectrum, lo, hi);
944
945	TF1 f("f", "[1]*(x/500)^[0]", lo, hi);
946	f.SetParameter(0, -3.0);
947	f.SetParameter(1, spectrum.GetMaximum());
948	f.SetLineColor(kBlue);
949	f.SetLineWidth(2);
950	spectrum.Fit(&f, "NIR"); // M skipped
951	f.DrawCopy("same");
952
953	TString str = FormString(f);
954
955	TLatex tex;
956	tex.SetTextSize(0.045);
957	tex.SetBit(TLatex::kTextNDC);
958	tex.SetTextAlign(31);
959	tex.DrawLatex(0.89, 0.935, str);
960
961	str = FormString(f, 1);
962	tex.DrawLatex(0.89, 0.83, str);
963
964	str = FormString(f, 2);
965	tex.DrawLatex(0.89, 0.735, str);
966
967	TArrayD res(2);
968	res[0] = f.GetParameter(0);
969	res[1] = f.GetParameter(1);
970	return res;
971	}
972
973	// --------------------------------------------------------------------------
974	//
975	// Calculate the final spectrum from:
976	// - collection area
977	// - excess
978	// - correction coefficients
979	// - ontime
980	// and display it
981	//
982	TArrayD MJSpectrum::DisplaySpectrum(MHCollectionArea &area, TH1D &excess, MHEnergyEst &hest, Double_t ontime) const
983	{
984	// Create copies of the histograms
985	TH1D collarea(area.GetHEnergy());
986	TH1D spectrum(excess);
987	TH1D weights;
988
989	// Get spill-over corrections from energy estimation
990	hest.GetWeights(weights); // E_mc/E_est
991
992	// Print effective on-time
993	cout << "Effective On time: " << ontime << "s" << endl;
994
995	// Setup spectrum plot
996	spectrum.SetNameTitle("Preliminary", "N/sm^{2} versus Energy (before unfolding)");
997	spectrum.SetYTitle("N/sm^{2}");
998	spectrum.SetDirectory(NULL);
999	spectrum.SetBit(kCanDelete);
1000
1001	// Divide by collection are and on-time
1002	spectrum.Scale(1./ontime);
1003	spectrum.Divide(&collarea);
1004
1005	// Draw spectrum before applying spill-over corrections
1006	TCanvas &c1 = fDisplay->AddTab("Spectrum");
1007	c1.Divide(2,2);
1008	c1.cd(1);
1009	gPad->SetBorderMode(0);
1010	gPad->SetLogx();
1011	gPad->SetLogy();
1012	gPad->SetGridx();
1013	gPad->SetGridy();
1014	collarea.DrawCopy();
1015
1016	c1.cd(2);
1017	gPad->SetBorderMode(0);
1018	gPad->SetLogx();
1019	gPad->SetLogy();
1020	gPad->SetGridx();
1021	gPad->SetGridy();
1022	TH1D spec=(TH1D)spectrum.DrawCopy();
1023	//FitSpectrum(*spec);
1024
1025	c1.cd(3);
1026	gPad->SetBorderMode(0);
1027	gPad->SetLogx();
1028	gPad->SetLogy();
1029	gPad->SetGridx();
1030	gPad->SetGridy();
1031	weights.DrawCopy();
1032
1033	// Apply spill-over correction (done't take the errors into account)
1034	// They are supposed to be identical with the errors of the
1035	// collection area and cancel out.
1036	//spectrum.Multiply(&weights);
1037	spectrum.SetNameTitle("Flux", "Spectrum");
1038	spectrum.SetBit(TH1::kNoStats);
1039
1040	// Minimum number of excessevents to get 3sigma in 1h
1041	TF1 sensl("SensLZA", "85*(x/200)^(-0.55)", 100, 6000);
1042	TF1 sensh("SensHZA", "35*(x/200)^(-0.70)", 100, 1000);
1043	//sens.SetLineColor(kBlue);
1044	//sens.DrawClone("Csame");
1045
1046	c1.cd(4);
1047	gPad->SetBorderMode(0);
1048	gPad->SetLogx();
1049	gPad->SetLogy();
1050	gPad->SetGridx();
1051	gPad->SetGridy();
1052
1053	TGraph gsensl;//("Sensitivity LZA", "");
1054	TGraph gsensh;//("Sensitivity HZA", "");
1055
1056	const Double_t sqrton = TMath::Sqrt(ontime/3600.);
1057
1058	for (int i=0; i<excess.GetNbinsX(); i++)
1059	{
1060	spectrum.SetBinContent(i+1, spectrum.GetBinContent(i+1)*weights.GetBinContent(i+1));
1061	spectrum.SetBinError(i+1, spectrum.GetBinError(i+1) *weights.GetBinContent(i+1));
1062
1063	spectrum.SetBinContent(i+1, spectrum.GetBinContent(i+1)/spectrum.GetBinWidth(i+1)*1000);
1064	spectrum.SetBinError(i+1, spectrum.GetBinError(i+1)/ spectrum.GetBinWidth(i+1)*1000);
1065
1066	if (collarea.GetBinContent(i+1)<=0)
1067	continue;
1068
1069	const Double_t cen = spectrum.GetBinCenter(i+1);
1070	gsensl.SetPoint(gsensl.GetN(), cen, sensl.Eval(cen)sqrton/spectrum.GetBinWidth(i+1)1000/collarea.GetBinContent(i+1)/ontime);
1071	gsensh.SetPoint(gsensh.GetN(), cen, sensh.Eval(cen)sqrton/spectrum.GetBinWidth(i+1)1000/collarea.GetBinContent(i+1)/ontime);
1072
1073	cout << cen << " " << sensl.Eval(cen)sqrton/spectrum.GetBinWidth(i+1)1000/collarea.GetBinContent(i+1)/ontime;
1074	cout << " " << sensh.Eval(cen)sqrton/spectrum.GetBinWidth(i+1)1000/collarea.GetBinContent(i+1)/ontime;
1075	cout << endl;
1076	}
1077
1078	spectrum.SetMinimum(1e-12);
1079	spectrum.SetXTitle("E [GeV]");
1080	spectrum.SetYTitle("dN/dE [N/TeVsm^{2}]");
1081	spec = (TH1D*)spectrum.DrawCopy();
1082
1083	TLatex tex;
1084	tex.SetTextColor(13);
1085
1086	TF1 fc("Crab", "6.0e-6(x/300)^(-2.31-0.26log10(x/300))", 100, 6000);
1087	fc.SetLineStyle(9);
1088	fc.SetLineWidth(2);
1089	fc.SetLineColor(14);
1090	fc.DrawCopy("same");
1091
1092	tex.DrawLatex(1250, fc.Eval(1250), "Crab/\\Gamma=-2.3");
1093
1094	TF1 fa("PG1553", "1.8e-6*(x/200)^-4.21", 90, 600);
1095	static_cast<const TAttLine&>(fc).Copy(fa);
1096	fa.DrawCopy("same");
1097
1098	tex.SetTextAlign(23);
1099	tex.DrawLatex(600, fa.Eval(600), "PG1553/\\Gamma=-4.2");
1100
1101	gROOT->SetSelectedPad(0);
1102
1103	gsensl.SetLineStyle(5);
1104	gsensl.SetLineColor(14);
1105	gsensl.SetLineWidth(2);
1106	gsensl.DrawClone("C")->SetBit(kCanDelete);
1107
1108	gsensh.SetLineStyle(5);
1109	gsensh.SetLineColor(14);
1110	gsensh.SetLineWidth(2);
1111	gsensh.DrawClone("C")->SetBit(kCanDelete);
1112
1113	// Display dN/dE*E^2 for conveinience
1114	fDisplay->AddTab("Check");
1115	gPad->SetLogx();
1116	gPad->SetLogy();
1117	gPad->SetGrid();
1118
1119	// Calculate Spectrum * E^2
1120	for (int i=0; i<spectrum.GetNbinsX(); i++)
1121	{
1122	const Double_t e = TMath::Sqrt(spectrum.GetBinLowEdge(i+1)spectrum.GetBinLowEdge(i+2))1e-3;
1123
1124	spectrum.SetBinContent(i+1, spectrum.GetBinContent(i+1)ee);
1125	spectrum.SetBinError( i+1, spectrum.GetBinError(i+1) ee);
1126	}
1127
1128	for (int i=0; i<gsensl.GetN(); i++)
1129	{
1130	const Double_t e = gsensl.GetX()[i]*1e-3;
1131
1132	gsensl.GetY()[i] = ee;
1133	gsensh.GetY()[i] = ee;
1134	}
1135
1136	spectrum.SetName("FluxStd");
1137	spectrum.SetMarkerStyle(kFullDotMedium);
1138	spectrum.SetTitle("Differential flux times E^{2}");
1139	spectrum.SetYTitle("E^{2}·dN/dE [N·TeV/sm^{2}]");
1140	spectrum.SetDirectory(0);
1141	spectrum.DrawCopy();
1142
1143	TF1 fc2("CrabE^2", "x^2Crab*1e-6", 100, 6000);
1144	static_cast<const TAttLine&>(fc).Copy(fc2);
1145	fc2.DrawCopy("same");
1146
1147	TF1 fa2("PG1553E^2", "x^2PG1553*1e-6", 100, 6000);
1148	static_cast<const TAttLine&>(fc).Copy(fa2);
1149	fa2.DrawCopy("same");
1150
1151	gsensl.DrawClone("C")->SetBit(kCanDelete);
1152	gsensh.DrawClone("C")->SetBit(kCanDelete);
1153
1154	// Fit Spectrum
1155	c1.cd(4);
1156	return FitSpectrum(*spec);
1157
1158	/*
1159	TF1 f("f", "[1]*(x/1e3)^[0]", 10, 3e4);
1160	f.SetParameter(0, -2.87);
1161	f.SetParameter(1, 1.9e-6);
1162	f.SetLineColor(kGreen);
1163	spectrum.Fit(&f, "NIM", "", 100, 5000);
1164	f.DrawCopy("same");
1165
1166	const Double_t p0 = f.GetParameter(0);
1167	const Double_t p1 = f.GetParameter(1);
1168
1169	const Double_t e0 = f.GetParError(0);
1170	const Double_t e1 = f.GetParError(1);
1171
1172	const Int_t np = TMath::Nint(TMath::Floor(TMath::Log10(p1)));
1173	const Double_t exp = TMath::Power(10, np);
1174
1175	TString str;
1176	str += Form("(%.2f#pm%.2f)10^{%d}", p1/exp, e1/exp, np);
1177	str += Form("(\\frac{E}{TeV})^{%.2f#pm%.2f}", p0, e0);
1178	str += "\\frac{ph}{TeVm^{2}s}";
1179
1180	TLatex tex;
1181	tex.SetTextSize(0.045);
1182	tex.SetBit(TLatex::kTextNDC);
1183	tex.DrawLatex(0.45, 0.935, str);
1184
1185	str = Form("\\chi^{2}/NDF=%.2f", f.GetChisquare()/f.GetNDF());
1186	tex.DrawLatex(0.70, 0.83, str);
1187
1188	TArrayD res(2);
1189	res[0] = f.GetParameter(0);
1190	res[1] = f.GetParameter(1);
1191
1192	return res;
1193	*/
1194	}
1195
1196	// --------------------------------------------------------------------------
1197	//
1198	// Scale some image parameter plots using the scale factor and plot them
1199	// together with the corresponding MC histograms.
1200	// Called from DisplaySize
1201	//
1202	Bool_t MJSpectrum::PlotSame(MStatusArray &arr, MParList &plist, const char name, const char tab, const char *plot, Double_t scale) const
1203	{
1204	TString same(name);
1205	same += "Same";
1206
1207	TH1 h1 = (TH1)arr.FindObjectInCanvas(name, "TH1F", tab);
1208	TH1 h2 = (TH1)arr.FindObjectInCanvas(same, "TH1F", tab);
1209	if (!h1 \|\| !h2)
1210	return kFALSE;
1211
1212	TObject *obj = plist.FindObject(plot);
1213	if (!obj)
1214	{
1215	*fLog << warn << plot << " not in parameter list... skipping." << endl;
1216	return kFALSE;
1217	}
1218
1219	TH1 h3 = (TH1)obj->FindObject(name);
1220	if (!h3)
1221	{
1222	*fLog << warn << name << " not found in " << plot << "... skipping." << endl;
1223	return kFALSE;
1224	}
1225
1226
1227	const MAlphaFitter fit = (MAlphaFitter)plist.FindObject("MAlphaFitter");
1228	const Double_t ascale = fit ? fit->GetScaleFactor() : 1;
1229
1230	gPad->SetBorderMode(0);
1231	h2->SetLineColor(kBlack);
1232	h3->SetLineColor(kBlue);
1233	h2->Add(h1, -ascale);
1234
1235	//h2->Scale(1./ontime); //h2->Integral());
1236	h3->Scale(scale); //h3->Integral());
1237
1238	h2->SetMaximum(1.05*TMath::Max(h2->GetMaximum(), h3->GetMaximum()));
1239
1240	h2 = h2->DrawCopy();
1241	h3 = h3->DrawCopy("same");
1242
1243	// Don't do this on the original object!
1244	h2->SetStats(kFALSE);
1245	h3->SetStats(kFALSE);
1246
1247	return kTRUE;
1248	}
1249
1250	// --------------------------------------------------------------------------
1251	//
1252	// Take a lot of histograms and plot them together in one plot.
1253	// Calls PlotSame
1254	//
1255	Bool_t MJSpectrum::DisplaySize(MParList &plist, Double_t scale) const
1256	{
1257	*fLog << inf << "Reading from file: " << fPathIn << endl;
1258
1259	TFile file(fPathIn, "READ");
1260	if (!file.IsOpen())
1261	{
1262	*fLog << err << dbginf << "ERROR - Could not open file " << fPathIn << endl;
1263	return kFALSE;
1264	}
1265
1266	file.cd();
1267	MStatusArray arr;
1268	if (arr.Read()<=0)
1269	{
1270	*fLog << "MStatusDisplay not found in file... abort." << endl;
1271	return kFALSE;
1272	}
1273
1274	TH1 excess = (TH1D)arr.FindObjectInCanvas("Excess", "TH1D", "Hist");
1275	if (!excess)
1276	return kFALSE;
1277
1278	// ------------------- Plot excess versus size -------------------
1279
1280	TCanvas &c = fDisplay->AddTab("Excess");
1281	c.Divide(3,2);
1282	c.cd(1);
1283	gPad->SetBorderMode(0);
1284	gPad->SetLogx();
1285	gPad->SetLogy();
1286	gPad->SetGridx();
1287	gPad->SetGridy();
1288
1289	excess->SetTitle("Excess events vs Size (data/black, mc/blue)");
1290	excess = excess->DrawCopy("E2");
1291	// Don't do this on the original object!
1292	excess->SetStats(kFALSE);
1293	excess->SetMarkerStyle(kFullDotMedium);
1294	excess->SetFillColor(kBlack);
1295	excess->SetFillStyle(0);
1296	excess->SetName("Excess ");
1297	excess->SetDirectory(0);
1298
1299	TObject *o=0;
1300	if ((o=plist.FindObject("ExcessMC")))
1301	{
1302	TH1 histsel = (TH1F)o->FindObject("");
1303	if (histsel)
1304	{
1305	if (scale<0)
1306	scale = excess->Integral()/histsel->Integral();
1307
1308	histsel->Scale(scale);
1309	histsel->SetLineColor(kBlue);
1310	histsel->SetBit(kCanDelete);
1311	histsel = histsel->DrawCopy("E1 same");
1312	// Don't do this on the original object!
1313	histsel->SetStats(kFALSE);
1314	/*
1315	fLog->Separator("Kolmogorov Test");
1316	histsel->KolmogorovTest(excess, "DX");
1317	fLog->Separator("Chi^2 Test");
1318	const Double_t p = histsel->Chi2Test(excess, "P");
1319
1320	TLatex tex;
1321	tex.SetBit(TLatex::kTextNDC);
1322	tex.DrawLatex(0.75, 0.93, Form("P(\\chi^{2})=%.0f%%", p*100));
1323	*/
1324	}
1325	}
1326
1327	// -------------- Comparison of Image Parameters --------------
1328	c.cd(2);
1329	PlotSame(arr, plist, "Dist", "HilSrc", "MHHilSrcMCPost", scale);
1330
1331	c.cd(3);
1332	PlotSame(arr, plist, "Length", "PostCut", "MHHillasMCPost", scale);
1333
1334	c.cd(4);
1335	PlotSame(arr, plist, "M3l", "HilExt", "MHHilExtMCPost", scale);
1336
1337	c.cd(5);
1338	PlotSame(arr, plist, "Conc1", "NewPar", "MHNewParMCPost", scale);
1339
1340	c.cd(6);
1341	PlotSame(arr, plist, "Width", "PostCut", "MHHillasMCPost", scale);
1342
1343	return kTRUE;
1344	}
1345
1346	// --------------------------------------------------------------------------
1347	//
1348	void MJSpectrum::DisplayCutEfficiency(const MHCollectionArea &area0, const MHCollectionArea &area1) const
1349	{
1350	if (!fDisplay)
1351	return;
1352
1353	const TH1D &trig = area0.GetHEnergy();
1354	TH1D &cut = (TH1D&)*area1.GetHEnergy().Clone();
1355
1356	fDisplay->AddTab("CutEff");
1357
1358	gPad->SetBorderMode(0);
1359	gPad->SetFrameBorderMode(0);
1360	gPad->SetLogx();
1361	gPad->SetGridx();
1362	gPad->SetGridy();
1363
1364	cut.Divide(&trig);
1365	cut.Scale(100);
1366	cut.SetNameTitle("CutEff", "Background Supression: Cut efficiency (after star)");
1367	cut.SetYTitle("\\eta [%]");
1368	cut.SetDirectory(0);
1369	cut.SetMinimum(0);
1370	cut.SetMaximum(100);
1371	cut.SetBit(kCanDelete);
1372	cut.Draw();
1373
1374	TLine line;
1375	line.SetLineColor(kBlue);
1376	line.SetLineWidth(2);
1377	line.SetLineStyle(kDashed);
1378	line.DrawLine(cut.GetBinLowEdge(1), 50, cut.GetBinLowEdge(cut.GetNbinsX()+1), 50);
1379	}
1380
1381	void MJSpectrum::SetupHistEvtDist(MHn &hist) const
1382	{
1383	hist.AddHist("MMcEvt.fEnergy");
1384	hist.InitName("EnergyDist;EnergyEst");
1385	hist.InitTitle("Unweighted event distribution (Real statistics);E [GeV];Counts;");
1386
1387	hist.AddHist("MPointingPos.fZd");
1388	hist.InitName("ZdDist;Theta");
1389	hist.InitTitle("Unweighted event distribution (Real statistics);Zd [\\circ];Counts;");
1390	}
1391
1392	void MJSpectrum::SetupHistEnergyEst(MHn &hist) const
1393	{
1394	const char *res = "log10(MEnergyEst.fVal)-log10(MMcEvt.fEnergy)";
1395
1396	hist.AddHist("MHillas.fSize", res);
1397	hist.InitName("ResSize;Size;EnergyResidual");
1398	hist.InitTitle(";S [phe];\\Delta lg E;");
1399	hist.SetDrawOption("colz profx");
1400
1401	hist.AddHist("MPointingPos.fZd", res);
1402	hist.InitName("ResTheta;Theta;EnergyResidual");
1403	hist.InitTitle(";Zd [\\circ];\\Delta lg E;");
1404	hist.SetDrawOption("colz profx");
1405
1406	hist.AddHist("MNewImagePar.fLeakage1", res);
1407	hist.InitName("ResLeak;Leakage;EnergyResidual");
1408	hist.InitTitle(";Leak;\\Delta lg E;");
1409	hist.SetDrawOption("colz profx");
1410
1411	hist.AddHist("MHillasSrc.fDist*3.37e-3", res);
1412	hist.InitName("ResDist;Dist;EnergyResidual");
1413	hist.InitTitle(";D [\\circ];\\Delta lg E;");
1414	hist.SetDrawOption("colz profx");
1415	}
1416
1417	void MJSpectrum::SetupHistDisp(MHn &hist) const
1418	{
1419	const char res = "-Disp.fValsign(MHillasSrc.fCosDeltaAlpha)-MHillasSrc.fDist*3.37e-3";
1420
1421	hist.AddHist("MHillas.fSize", res);
1422	hist.InitName("ResSize;Size;ResidualDist");
1423	hist.InitTitle(";S [phe];Disp-Dist [\\circ];");
1424	hist.SetDrawOption("colz profx");
1425
1426	hist.AddHist("MPointingPos.fZd", res);
1427	hist.InitName("ResTheta;Theta;ResidualDist");
1428	hist.InitTitle(";Zd [\\circ];Disp-Dist [\\circ];");
1429	hist.SetDrawOption("colz profx");
1430
1431	hist.AddHist("MNewImagePar.fLeakage1", res);
1432	hist.InitName("ResLeak;Leakage;ResidualDist");
1433	hist.InitTitle(";Leak;Disp-Dist [\\circ];");
1434	hist.SetDrawOption("colz profx");
1435
1436	hist.AddHist("MHillasExt.fSlopeLong*sign(MHillasSrc.fCosDeltaAlpha)/3.37e-3", res);
1437	hist.InitName("ResSlope;Slope;ResidualDist");
1438	hist.InitTitle(";Slope;Disp-Dist [\\circ];");
1439	hist.SetDrawOption("colz profx");
1440	}
1441
1442	void MJSpectrum::SetupHistEnergyRes(MHn &hist) const
1443	{
1444	hist.AddHist("MEnergyEst.fVal", "(MMcEvt.fEnergy/MEnergyEst.fVal-1)^2", MH3::kProfile);
1445	hist.InitName("ResEest;EnergyEst;");
1446	hist.InitTitle(";E_{est} [GeV];Resolution (E_{mc}/E_{est}-1)^{2};");
1447
1448	//hist.AddHist("MMcEvt.fEnergy", "(MEnergyEst.fVal/MMcEvt.fEnergy-1)^2", MH3::kProfile);
1449	//hist.InitName("ResEmc;EnergyEst;");
1450	//hist.InitTitle(";E_{mc} [GeV];Resolution (E_{est}/E_{mc}-1)^{2};");
1451	hist.AddHist("MHillas.fSize", "(MMcEvt.fEnergy/MEnergyEst.fVal-1)^2", MH3::kProfile);
1452	hist.InitName("ResSize;Size;");
1453	hist.InitTitle(";S [phe];Resolution (E_{mc}/E_{est}-1)^{2};");
1454
1455	hist.AddHist("MPointingPos.fZd", "(MMcEvt.fEnergy/MEnergyEst.fVal-1)^2", MH3::kProfile);
1456	hist.InitName("ResTheta;Theta;");
1457	hist.InitTitle(";\\Theta [\\circ];Resolution (E_{mc}/E_{est}-1)^{2};");
1458
1459	hist.AddHist("MMcEvt.fImpact/100", "(MMcEvt.fEnergy/MEnergyEst.fVal-1)^2", MH3::kProfile);
1460	hist.InitName("ResImpact;Impact;");
1461	hist.InitTitle(";I [m];Resolution (E_{mc}/E_{est}-1)^{2};");
1462	}
1463
1464	// --------------------------------------------------------------------------
1465	//
1466	// Setup write to write:
1467	// container tree optional?
1468	// -------------- ---------- -----------
1469	// "MHillas" to "Events"
1470	// "MHillasSrc" to "Events"
1471	// "Hadronness" to "Events" yes
1472	// "MEnergyEst" to "Events" yes
1473	// "DataType" to "Events"
1474	//
1475	void MJSpectrum::SetupWriter(MWriteRootFile write/, const char name/) const
1476	{
1477	if (!write)
1478	return;
1479
1480	//write->SetName(name);
1481	write->AddContainer("MHillas", "Events");
1482	write->AddContainer("MHillasSrc", "Events");
1483	write->AddContainer("MHillasExt", "Events");
1484	//write->AddContainer("MPointingPos", "Events");
1485	write->AddContainer("MHillasSrcAnti", "Events", kFALSE);
1486	write->AddContainer("MImagePar", "Events", kFALSE);
1487	write->AddContainer("MNewImagePar", "Events", kFALSE);
1488	write->AddContainer("MNewImagePar2", "Events", kFALSE);
1489	write->AddContainer("Hadronness", "Events", kFALSE);
1490	write->AddContainer("MSrcPosCam", "Events", kFALSE);
1491	write->AddContainer("MSrcPosAnti", "Events", kFALSE);
1492	write->AddContainer("ThetaSquared", "Events", kFALSE);
1493	write->AddContainer("OpticalAxis", "Events", kFALSE);
1494	write->AddContainer("Disp", "Events", kFALSE);
1495	write->AddContainer("Ghostbuster", "Events", kFALSE);
1496	write->AddContainer("MEnergyEst", "Events", kFALSE);
1497	write->AddContainer("MTime", "Events", kFALSE);
1498	write->AddContainer("MMcEvt", "Events", kFALSE);
1499	write->AddContainer("MWeight", "Events");
1500	write->AddContainer("DataType", "Events");
1501	write->AddContainer("FileId", "Events");
1502	write->AddContainer("EvtNumber", "Events");
1503	}
1504
1505	Bool_t MJSpectrum::Process(const MDataSet &set)
1506	{
1507	if (!set.IsValid())
1508	{
1509	*fLog << err << "ERROR - DataSet invalid!" << endl;
1510	return kFALSE;
1511	}
1512
1513	if (!HasWritePermission(GetPathOut()))
1514	return kFALSE;
1515
1516	if (!CheckEnv())
1517	return kFALSE;
1518
1519	// --------------------------------------------------------------------------------
1520
1521	*fLog << inf;
1522	fLog->Separator(GetDescriptor());
1523	*fLog << "Compile Monte Carlo sample (dataset " << set.GetBaseName() << ")" << endl;
1524	*fLog << endl;
1525
1526	if (fDisplay)
1527	fDisplay->SetWindowName(fName);
1528
1529	// Setup everything which is read from the ganymed file
1530	MBinning bins1("BinningAlpha");
1531	MBinning bins2("BinningEnergyEst");
1532	MBinning bins3("BinningTheta");
1533	MBinning bins4("BinningFalseSource");
1534	MBinning bins5("BinningWidth");
1535	MBinning bins6("BinningLength");
1536	MBinning bins7("BinningDist");
1537	MBinning bins8("BinningM3Long");
1538	MBinning bins9("BinningM3Trans");
1539	MBinning bins0("BinningSlope");
1540	MBinning binsa("BinningAsym");
1541	MBinning binsb("BinningConc1");
1542
1543	MEnv env("", "ganymed.rc");
1544
1545	MAlphaFitter fit;
1546
1547	MParList plist;
1548	plist.AddToList(&bins0);
1549	plist.AddToList(&bins1);
1550	plist.AddToList(&bins3);
1551	plist.AddToList(&bins4);
1552	plist.AddToList(&bins5);
1553	plist.AddToList(&bins6);
1554	plist.AddToList(&bins7);
1555	plist.AddToList(&bins8);
1556	plist.AddToList(&bins9);
1557	plist.AddToList(&binsa);
1558	plist.AddToList(&binsb);
1559	plist.AddToList(&fit);
1560
1561	// Read from the ganymed file
1562	TH1D htheta, size;
1563	Float_t ontime = ReadInput(plist, htheta, size);
1564	if (ontime<0)
1565	{
1566	*fLog << err << GetDescriptor() << ": Could not determine effective on time..." << endl;
1567	return kFALSE;
1568	}
1569
1570	// Set Zenith angle binning to binning from the ganymed-file
1571	bins3.SetEdges(htheta, 'x');
1572
1573	// Read energy binning from resource file
1574	if (!CheckEnv(bins2))
1575	{
1576	*fLog << err << "ERROR - Reading energy binning from resources failed." << endl;
1577	return kFALSE;
1578	}
1579	plist.AddToList(&bins2); // For later use in MC processing
1580
1581	// -------------- Fill excess events versus energy ---------------
1582
1583	TH1D excess;
1584
1585	if (fForceOnTimeFit)
1586	{
1587	// Refill excess histogram to determine the excess events
1588	// If we use the eff. on-time fit we have to loop over the data first
1589	// This is not really desired, because if something is wrong with
1590	// the Monte Carlos the program runs quite long before it fails
1591	if (!Refill(plist, excess))
1592	return kFALSE;
1593
1594	// Print the setup and result of the MAlphaFitter, print used cuts
1595	PrintSetup(fit);
1596
1597	// ------------ On user request redo eff. on-time fit ------------
1598	const MHEffectiveOnTime htime = (MHEffectiveOnTime)plist.FindObject("MHEffectiveOnTime");
1599	if (!htime)
1600	{
1601	// This should never happen, bt you never know
1602	*fLog << err;
1603	*fLog << "ERROR - Use of new effective on-time fit requested for on-time determination," << endl;
1604	*fLog << " but MHEffectiveOnTime not found in parameter list... aborting." << endl;
1605	return kFALSE;
1606	}
1607
1608	const TH1D &h = htime->GetHEffOnTheta();
1609	/*
1610	if (!htime->IsConsistent() \|\| h.GetNbinsX()!=htheta.GetNbinsX())
1611	{
1612	*fLog << err << "ERROR - Effective on-time from newly filles MHEffectiveOnTime (Tab='OnTime) invalid... aborting." << endl;
1613	return kFALSE;
1614	}*/
1615
1616	*fLog << inf;
1617	*fLog << "Using eff. on-time from new MHEffectiveOnTime (see also 'OnTime')" << endl;
1618	*fLog << " Orig. value: " << ontime << "s" << endl;
1619
1620	// Copy ontime from newly filled and fitted eff on-time histogram
1621	ontime = htime->GetEffOnTime();
1622
1623	*fLog << " New value: " << ontime << "s" << endl;
1624
1625	h.Copy(htheta); // Copy contents of newly filled hist into on-time vs. theta
1626	htheta.SetName("Theta"); // Copy overwrites the name needed in DisplayResult
1627	htheta.SetDirectory(0); // Remove from global directory added by SetName
1628	}
1629
1630	// ---------------------------------------------------------------
1631
1632	// Initialize weighting to a new spectrum as defined in the resource file
1633	MMcSpectrumWeight weight;
1634	if (!InitWeighting(set, weight))
1635	return kFALSE;
1636
1637	// Print Theta and energy binning for cross-checks
1638	*fLog << all << endl;
1639	bins2.Print();
1640	bins3.Print();
1641
1642	// Now we read the MC distribution as produced by corsika
1643	// vs zenith angle and energy.
1644	// Weight for the new energy spectrum defined in MMcSpectumWeight
1645	// are applied.
1646	// Also correction for different lower energy bounds and
1647	// different production areas (impact parameters) are applied.
1648	TH2D hist;
1649	hist.UseCurrentStyle();
1650	MH::SetBinning(hist, bins3, bins2);
1651	if (!ReadOrigMCDistribution(set, hist, weight))
1652	return kFALSE;
1653
1654	// Check if user has closed the display
1655	if (!fDisplay)
1656	return kTRUE;
1657
1658	// Display histograms and calculate za-weights into htheta
1659	GetThetaDistribution(htheta, hist);
1660
1661	// Give the zenith angle weights to the weighting task
1662	weight.SetWeightsZd(&htheta);
1663
1664	// No we apply the the zenith-angle-weights to the corsika produced
1665	// MC distribution. Unfortunately this must be done manually
1666	// because we are multiplying column by column
1667	for (int y=0; y<=hist.GetNbinsY()+1; y++)
1668	for (int x=0; x<=hist.GetNbinsX()+1; x++)
1669	{
1670	hist.SetBinContent(x, y, hist.GetBinContent(x, y)*htheta.GetBinContent(x));
1671	hist.SetBinError(x, y, hist.GetBinError(x, y) *htheta.GetBinContent(x));
1672	}
1673
1674	// Display the resulting distribution and check it matches
1675	// the observation time distribution (this could be false
1676	// for example if you miss MCs of some zenith angles, which you have
1677	// data for)
1678	if (!DisplayResult(hist))
1679	return kFALSE;
1680
1681	// Refill excess histogram to determine the excess events
1682	if (!fForceOnTimeFit)
1683	{
1684	if (!Refill(plist, excess))
1685	return kFALSE;
1686
1687	// Print the setup and result of the MAlphaFitter, print used cuts
1688	PrintSetup(fit);
1689	}
1690
1691	// ------------------------- Final loop --------------------------
1692
1693	*fLog << endl;
1694	fLog->Separator("Calculate efficiencies");
1695	*fLog << endl;
1696
1697	MTaskList tlist2;
1698	plist.AddToList(&tlist2);
1699
1700	MReadMarsFile read("Events");
1701	read.DisableAutoScheme();
1702	if (!set.AddFilesOn(read))
1703	return kFALSE;
1704
1705	// Selector to get correct (final) theta-distribution
1706	//temp1.SetXTitle("MPointingPos.fZd");
1707	//
1708	//MH3 mh3(temp1);
1709	//
1710	//MFEventSelector2 sel2(mh3);
1711	//sel2.SetHistIsProbability();
1712	//
1713	//MContinue contsel(&sel2);
1714	//contsel.SetInverted();
1715
1716	// Get correct source position
1717	//MSrcPosCalc calc;
1718
1719	// Calculate corresponding Hillas parameters
1720	//MHillasCalc hcalc1;
1721	//MHillasCalc hcalc2("MHillasCalcAnti");
1722	//hcalc1.SetFlags(MHillasCalc::kCalcHillasSrc);
1723	//hcalc2.SetFlags(MHillasCalc::kCalcHillasSrc);
1724	//hcalc2.SetNameHillasSrc("MHillasSrcAnti");
1725	//hcalc2.SetNameSrcPosCam("MSrcPosAnti");
1726
1727	// Fill collection area and energy estimator (unfolding)
1728	// Make sure to use the same binning for MHCollectionArea and MHEnergyEst
1729	MHCollectionArea area0("TriggerArea");
1730	MHCollectionArea area1;
1731	area0.SetHistAll(hist);
1732	area1.SetHistAll(hist);
1733
1734	MHEnergyEst hest;
1735
1736	MFillH fill30(&area0, "", "FillTriggerArea");
1737	MFillH fill31(&area1, "", "FillCollectionArea");
1738	MFillH fill4(&hest, "", "FillEnergyEst");
1739	MFillH fill5("MHThreshold", "", "FillThreshold");
1740	fill30.SetWeight();
1741	fill31.SetWeight();
1742	fill4.SetWeight();
1743	fill5.SetWeight();
1744	fill30.SetNameTab("TrigArea");
1745	fill31.SetNameTab("ColArea");
1746	fill4.SetNameTab("E-Est");
1747	fill5.SetNameTab("Threshold");
1748
1749	/*
1750	MH3 henergy("MMcEvt.fEnergy");
1751	henergy.SetName("EventDist;EnergyEst");
1752	henergy.SetTitle("Unweighted event distribution (Real statistics);E [GeV];Counts;");
1753	henergy.Sumw2();
1754	*/
1755
1756	// ---------------------------------------------------------
1757
1758	MBinning binsA(50, 10, 100000, "BinningSize", "log");
1759	MBinning binsC(50, 0, 0.3, "BinningLeakage", "lin");
1760	MBinning binsB(51, -1, 1, "BinningEnergyResidual", "lin");
1761	MBinning binsD(51, -1, 1, "BinningResidualDist", "lin");
1762	MBinning binsI(16, 0, 800, "BinningImpact", "lin");
1763
1764	plist.AddToList(&binsA);
1765	plist.AddToList(&binsB);
1766	plist.AddToList(&binsC);
1767	plist.AddToList(&binsD);
1768	plist.AddToList(&binsI);
1769
1770	MHn heest("Energy", "Energy Residual (lg E_{est} - lg E_{mc})");
1771	SetupHistEnergyEst(heest);
1772
1773	MHn hdisp("Disp", "Dist residual (Disp-Dist)");
1774	SetupHistDisp(hdisp);
1775
1776	MHn henergy("EvtDist");
1777	SetupHistEvtDist(henergy);
1778
1779	MHn heres("EnergyRes");
1780	SetupHistEnergyRes(heres);
1781
1782	MFillH fill4b(&heest, "", "FillEnergyResidual");
1783	fill4b.SetWeight();
1784
1785	MFillH fill4c(&hdisp, "", "FillDispResidual");
1786	fill4c.SetWeight();
1787
1788	MFillH fill4d(&heres, "", "FillEnergyResolution");
1789	fill4d.SetWeight();
1790
1791	MFDataPhrase fdisp("Disp.fVal*sign(MHillasSrc.fCosDeltaAlpha)<0", "FilterDisp");
1792	fill4c.SetFilter(&fdisp);
1793
1794	// ---------------------------------------------------------
1795
1796	MH3 hsize("MHillas.fSize");
1797	hsize.SetName("ExcessMC");
1798	hsize.Sumw2();
1799
1800	MBinning bins(size, "BinningExcessMC");
1801	plist.AddToList(&hsize);
1802	plist.AddToList(&bins);
1803
1804	// ---------------------------------------------------------
1805
1806	MFillH fillsp("MHSrcPosCam", "MSrcPosCam", "FillSrcPosCam");
1807	MFillH fill0a(&henergy, "", "FillEventDist");
1808	MFillH fill1a("MHHillasMCPre [MHHillas]", "MHillas", "FillHillasPre");
1809	MFillH fill2a("MHHillasMCPost [MHHillas]", "MHillas", "FillHillasPost");
1810	MFillH fill3a("MHVsSizeMCPost [MHVsSize]", "MHillasSrc", "FillVsSizePost");
1811	MFillH fill4a("MHHilExtMCPost [MHHillasExt]", "MHillasSrc", "FillHilExtPost");
1812	MFillH fill5a("MHHilSrcMCPost [MHHillasSrc]", "MHillasSrc", "FillHilSrcPost");
1813	MFillH fill6a("MHImgParMCPost [MHImagePar]", "MImagePar", "FillImgParPost");
1814	MFillH fill7a("MHNewParMCPost [MHNewImagePar]", "MNewImagePar", "FillNewParPost");
1815	MFillH fill8a("ExcessMC [MH3]", "", "FillExcessMC");
1816	fillsp.SetNameTab("SrcPos");
1817	fill0a.SetNameTab("EvtDist");
1818	fill1a.SetNameTab("PreCut");
1819	fill2a.SetNameTab("PostCut");
1820	fill3a.SetNameTab("VsSize");
1821	fill4a.SetNameTab("HilExt");
1822	fill5a.SetNameTab("HilSrc");
1823	fill6a.SetNameTab("ImgPar");
1824	fill7a.SetNameTab("NewPar");
1825	fill8a.SetBit(MFillH::kDoNotDisplay);
1826	fill1a.SetWeight();
1827	fill2a.SetWeight();
1828	fill3a.SetWeight();
1829	fill4a.SetWeight();
1830	fill5a.SetWeight();
1831	fill6a.SetWeight();
1832	fill7a.SetWeight();
1833	fill8a.SetWeight();
1834	fillsp.SetWeight();
1835
1836	// FIXME: To be done: A task checking the lower 1% after the lower
1837	// energy limit!
1838
1839	MTaskEnv taskenv0("CalcDisp");
1840	taskenv0.SetDefault(fCalcDisp);
1841
1842	MTaskEnv taskenv1("CalcHadronness");
1843	taskenv1.SetDefault(fCalcHadronness);
1844
1845	MEnergyEstimate est;
1846	MTaskEnv taskenv2("EstimateEnergy");
1847	taskenv2.SetDefault(fEstimateEnergy ? fEstimateEnergy : &est);
1848
1849	MWriteRootFile write(GetPathOut());
1850	SetupWriter(&write);
1851
1852	MParameterI par = (MParameterI)plist.FindCreateObj("MParameterI", "DataType");
1853	if (!par)
1854	return kFALSE;
1855	par->SetVal(2);
1856
1857	// Not really necessary but for sanity
1858	TObject *cam = plist.FindObject("MSrcPosCam");
1859	if (cam)
1860	cam->Clear();
1861
1862	tlist2.AddToList(&read);
1863	// If no weighting should be applied but the events should
1864	// be thrown away according to the theta distribution
1865	// it is enabled here
1866	//if (!fRawMc && fNoThetaWeights)
1867	// tlist2.AddToList(&contsel);
1868	//tlist2.AddToList(&calc);
1869	//tlist2.AddToList(&hcalc1);
1870	//tlist2.AddToList(&hcalc2);
1871	tlist2.AddToList(&weight);
1872	tlist2.AddToList(&fillsp);
1873	tlist2.AddToList(&fill1a);
1874	tlist2.AddToList(&fill30);
1875	tlist2.AddToList(fCutQ);
1876	tlist2.AddToList(fCut0);
1877	tlist2.AddToList(&taskenv0);
1878	tlist2.AddToList(&taskenv1);
1879	tlist2.AddToList(&fdisp);
1880	tlist2.AddToList(&fill4c);
1881	tlist2.AddToList(fCut1);
1882	tlist2.AddToList(fCutS);
1883	tlist2.AddToList(fCut2);
1884	tlist2.AddToList(fCut3);
1885	tlist2.AddToList(&taskenv2);
1886	if (!GetPathOut().IsNull())
1887	tlist2.AddToList(&write);
1888	tlist2.AddToList(&fill31);
1889	tlist2.AddToList(&fill4);
1890	tlist2.AddToList(&fill4b);
1891	tlist2.AddToList(&fill4d);
1892	tlist2.AddToList(&fill5);
1893	tlist2.AddToList(&fill0a);
1894	tlist2.AddToList(&fill2a);
1895	tlist2.AddToList(&fill3a);
1896	tlist2.AddToList(&fill4a);
1897	tlist2.AddToList(&fill5a);
1898	tlist2.AddToList(&fill6a);
1899	tlist2.AddToList(&fill7a);
1900	tlist2.AddToList(&fill8a);
1901	//tlist2.AddToList(&fill9a);
1902
1903	// by setting it here it is distributed to all consecutive tasks
1904	tlist2.SetAccelerator(MTask::kAccDontReset\|MTask::kAccDontTime);
1905
1906	MEvtLoop loop2(fName); // *** fName ***
1907	loop2.SetParList(&plist);
1908	loop2.SetDisplay(fDisplay);
1909	loop2.SetLogStream(fLog);
1910
1911	if (!SetupEnv(loop2))
1912	return kFALSE;
1913
1914	if (!loop2.Eventloop(fMaxEvents))
1915	{
1916	*fLog << err << GetDescriptor() << ": Processing of MC-data failed." << endl;
1917	return kFALSE;
1918	}
1919
1920	if (!loop2.GetDisplay())
1921	{
1922	*fLog << err << GetDescriptor() << ": Execution stopped by user." << endl;
1923	return kFALSE;
1924	}
1925
1926	fLog->Separator("Energy Estimator");
1927	if (plist.FindObject("EstimateEnergy"))
1928	plist.FindObject("EstimateEnergy")->Print();
1929
1930	fLog->Separator("Spectrum");
1931
1932	// -------------------------- Spectrum ----------------------------
1933
1934	// Calculate and display spectrum (N/TeVsm^2 at 1TeV)
1935	TArrayD res(DisplaySpectrum(area1, excess, hest, ontime));
1936
1937	// Spectrum fitted (convert res[1] from TeV to GeV)
1938	TF1 flx("flx", MString::Format("%e*pow(x/500, %f)", res[1]/500, res[0]).Data());
1939
1940	// Number of events this spectrum would produce per s and m^2
1941	Double_t n = flx.Integral(weight.GetEnergyMin(), weight.GetEnergyMax());
1942
1943	// scale with effective collection area to get the event rate (N/s)
1944	// scale with the effective observation time to absolute observed events
1945	n = area1.GetCollectionAreaAbs()ontime; // N
1946
1947	// Now calculate the scale factor from the number of events
1948	// produced and the number of events which should have been
1949	// observed with our telescope in the time ontime
1950	const Double_t scale = n/area1.GetEntries();
1951
1952	// Print normalization constant
1953	cout << "MC normalization factor: " << scale << endl;
1954
1955	// Display cut efficiency
1956	DisplayCutEfficiency(area0, area1);
1957
1958	// Overlay normalized plots
1959	DisplaySize(plist, scale);
1960
1961	// check if output should be written
1962	if (!fPathOut.IsNull())
1963	{
1964	TNamed ganame("ganymed.root", fPathIn.Data());
1965	TNamed cmdline("CommandLine", fCommandLine.Data());
1966
1967	// Write the output
1968	TObjArray cont;
1969	cont.Add(&env); // ganymed.rc
1970	cont.Add(const_cast<TEnv*>(GetEnv())); // sponde.rc
1971	cont.Add(const_cast<MDataSet*>(&set)); // Dataset
1972	cont.Add(plist.FindObject("MAlphaFitter"));
1973	cont.Add(&area0);
1974	cont.Add(&area1);
1975	cont.Add(&hest);
1976	cont.Add(&ganame);
1977	cont.Add(&cmdline);
1978
1979	if (fDisplay)
1980	cont.Add(fDisplay);
1981
1982	if (!WriteContainer(cont, "", GetPathOut().IsNull()?"RECREATE":"UPDATE"))
1983	{
1984	*fLog << err << GetDescriptor() << ": Writing result failed." << endl;
1985	return kFALSE;
1986	}
1987	}
1988
1989	*fLog << all << GetDescriptor() << ": Done." << endl;
1990	*fLog << endl << endl;
1991
1992	return kTRUE;
1993	}

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

Download in other formats: