Context Navigation

source: trunk/MagicSoft/Mars/mjobs/MJCut.cc@ 7389

Visit:

Last change on this file since 7389 was 7389, checked in by tbretz, 19 years ago
* empty log message *
File size: 24.9 KB

Line
1	/* ======================================================================== *\
2	!
3	! *
4	! * This file is part of MARS, the MAGIC Analysis and Reconstruction
5	! * Software. It is distributed to you in the hope that it can be a useful
6	! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
7	! * It is distributed WITHOUT ANY WARRANTY.
8	! *
9	! * Permission to use, copy, modify and distribute this software and its
10	! * documentation for any purpose is hereby granted without fee,
11	! * provided that the above copyright notice appear in all copies and
12	! * that both that copyright notice and this permission notice appear
13	! * in supporting documentation. It is provided "as is" without express
14	! * or implied warranty.
15	! *
16	!
17	!
18	! Author(s): Thomas Bretz, 1/2005 <mailto:tbretz@astro.uni-wuerzburg.de>
19	!
20	! Copyright: MAGIC Software Development, 2000-2005
21	!
22	!
23	\* ======================================================================== */
24
25	/////////////////////////////////////////////////////////////////////////////
26	//
27	// MJCut
28	//
29	// FIXME: Preparation for wobble mode missing
30	//
31	/////////////////////////////////////////////////////////////////////////////
32	#include "MJCut.h"
33
34	// Root
35	#include <TEnv.h>
36	#include <TFile.h>
37
38	// Environment
39	#include "MLog.h"
40	#include "MLogManip.h"
41
42	// Eventloop
43	#include "MParList.h"
44	#include "MTaskList.h"
45	#include "MEvtLoop.h"
46
47	// Display
48	#include "MStatusDisplay.h"
49
50	// Tasks
51	#include "MReadReports.h"
52	#include "MReadMarsFile.h"
53	#include "MPrint.h"
54	#include "MContinue.h"
55	#include "MTaskEnv.h"
56	#include "MPointingDevCalc.h"
57	#include "MSrcPosCalc.h"
58	#include "MSrcPosCorrect.h"
59	#include "MHillasCalc.h"
60	#include "MFillH.h"
61	#include "MWriteRootFile.h"
62
63	// Fit signal environment
64	#include "../mhflux/MAlphaFitter.h"
65	#include "../mhflux/MHAlpha.h"
66
67	// Containers
68	#include "MH3.h"
69	#include "MBinning.h"
70	#include "MDataSet.h"
71	#include "MParameters.h"
72	#include "MPointingPos.h"
73	#include "MObservatory.h"
74
75	ClassImp(MJCut);
76
77	using namespace std;
78
79	// --------------------------------------------------------------------------
80	//
81	// Default constructor. Set defaults for fStoreSummary, fStoreresult,
82	// fWriteOnly, fIsWobble and fFullDisplay to kFALSE and initialize
83	// /fEstimateEnergy and/ fCalcHadronness with NULL.
84	//
85	MJCut::MJCut(const char name, const char title)
86	: fStoreSummary(kFALSE), fStoreResult(kTRUE), fWriteOnly(kFALSE),
87	fIsWobble(kFALSE), fIsMonteCarlo(kFALSE), fFullDisplay(kTRUE),
88	fNameHist("MHThetaSq"), fCalcHadronness(0)
89	{
90	fName = name ? name : "MJCut";
91	fTitle = title ? title : "Standard program to perform g/h-separation cuts";
92	}
93
94	// --------------------------------------------------------------------------
95	//
96	// Destructor. Delete fEstimateEnergy and fCalcHadronness if != NULL
97	//
98	MJCut::~MJCut()
99	{
100	//if (fEstimateEnergy)
101	// delete fEstimateEnergy;
102	if (fCalcHadronness)
103	delete fCalcHadronness;
104	}
105
106	// --------------------------------------------------------------------------
107	//
108	// Set the name of the summary file (events after cut0)
109	// If you give a name the storage of this file is enabled implicitly.
110	// If you give no filename the storage is neither enabled nor disabled,
111	// but the storage file name is reset.
112	// If no filename is set the default filename is used.
113	// You can explicitly enable or disable the storage using EnableStoreOf*()
114	// The default argument is no filename.
115	//
116	void MJCut::SetNameSummaryFile(const char *name)
117	{
118	fNameSummary=name;
119	if (!fNameSummary.IsNull())
120	EnableStorageOfSummary();
121	}
122
123	// --------------------------------------------------------------------------
124	//
125	// Set the name of the summary file (events after cut3)
126	// If you give a name the storage of this file is enabled implicitly.
127	// If you give no filename the storage is neither enabled nor disabled,
128	// but the storage file name is reset.
129	// If no filename is set the default filename is used.
130	// You can explicitly enable or disable the storage using EnableStoreOf*()
131	// The default argument is no filename.
132	//
133	void MJCut::SetNameResultFile(const char *name)
134	{
135	fNameResult=name;
136	if (!fNameResult.IsNull())
137	EnableStorageOfResult();
138	}
139
140	// --------------------------------------------------------------------------
141	//
142	// Setup a task estimating the energy. The given task is cloned.
143	//
144	/*
145	void MJCut::SetEnergyEstimator(const MTask *task)
146	{
147	if (fEstimateEnergy)
148	delete fEstimateEnergy;
149	fEstimateEnergy = task ? (MTask*)task->Clone() : 0;
150	}
151	*/
152
153	// --------------------------------------------------------------------------
154	//
155	// Setup a task calculating the hadronness. The given task is cloned.
156	//
157	void MJCut::SetHadronnessCalculator(const MTask *task)
158	{
159	if (fCalcHadronness)
160	delete fCalcHadronness;
161	fCalcHadronness = task ? (MTask*)task->Clone() : 0;
162	}
163
164	// --------------------------------------------------------------------------
165	//
166	// return fOutputPath+"/ganymed%08d.root", num
167	//
168	TString MJCut::GetOutputFile(UInt_t num) const
169	{
170	TString p(fPathOut);
171	p += "/";
172	p += fNameOutput.IsNull() ? Form("ganymed%08d.root", num) : fNameOutput.Data();
173	gSystem->ExpandPathName(p);
174	return p;
175	}
176
177	/*
178	Bool_t MJCut::ReadTasks(const char fname, MTask &env1, MTask* &env2) const
179	{
180	// const TString fname = Form("%s/calib%08d.root", fPathIn.Data(), fSequence.GetSequence());
181
182	*fLog << inf << "Reading from file: " << fname << endl;
183
184	TFile file(fname, "READ");
185	if (!file.IsOpen())
186	{
187	*fLog << err << dbginf << "ERROR - Could not open file " << fname << endl;
188	return kFALSE;
189	}
190
191	TObject *o = file.Get("EstimateEnergy");
192	if (o && !o->InheritsFrom(MTask::Class()))
193	{
194	*fLog << err << dbginf << "ERROR - EstimateEnergy read from " << fname << " doesn't inherit from MTask!" << endl;
195	return kFALSE;
196	}
197	env1 = o ? (MTask*)o->Clone() : FNULL;
198
199	o = file.Get("CalcHadronness");
200	if (o && !o->InheritsFrom(MTask::Class()))
201	{
202	*fLog << err << dbginf << "ERROR - CalcHadronness read from " << fname << " doesn't inherit from MTask!" << endl;
203	return kFALSE;
204	}
205	env2 = o ? (MTask*)o->Clone() : NULL;
206
207	return kTRUE;
208	}
209	*/
210
211	// --------------------------------------------------------------------------
212	//
213	// Write the tasks in cont to the file corresponding to analysis number num,
214	// see GetOutputFile()
215	//
216	Bool_t MJCut::WriteTasks(UInt_t num, TObjArray &cont) const
217	{
218	if (fPathOut.IsNull())
219	{
220	*fLog << inf << "No output path specified via SetPathOut - no output written." << endl;
221	return kTRUE;
222	}
223
224	const TString oname(GetOutputFile(num));
225
226	*fLog << inf << "Writing to file: " << oname << endl;
227
228	TFile *file = 0;
229	if (fNameResult.IsNull() && fStoreResult)
230	{
231	file = (TFile*)gROOT->GetListOfFiles()->FindObject(oname);
232	if (file)
233	file->cd();
234	}
235	else
236	file = TFile::Open(oname, fOverwrite?"RECREATE":"NEW", "File created by MJCut", 9);
237
238	if (!file)
239	{
240	*fLog << err << "ERROR - Couldn't open file " << oname << " for writing..." << endl;
241	return kFALSE;
242	}
243
244	const Bool_t rc = WriteContainer(cont);
245
246	if (!(fNameResult.IsNull() && fStoreResult))
247	delete file;
248
249	return rc;
250	}
251
252	// --------------------------------------------------------------------------
253	//
254	// Write the result plots and other results to the file corresponding to
255	// analysis number num, see GetOutputFile()
256	//
257	Bool_t MJCut::WriteResult(const MParList &plist, UInt_t num) const
258	{
259	TObjArray arr;
260
261	// Save all MBinnings
262	TIter Next(plist);
263	TObject *o=0;
264	while ((o=Next()))
265	if (o->InheritsFrom(MBinning::Class()))
266	arr.Add(o);
267
268	// Save also the result, not only the setup
269	const MHAlpha halpha = (MHAlpha)plist.FindObject("Hist", "MHAlpha");
270	if (halpha)
271	arr.Add((TObject*)(&halpha->GetAlphaFitter()));
272
273	const TString fname(fNameOutput.IsNull() ? Form("ganymed%08d.root", num) : fNameOutput.Data());
274
275	// If requested, write to already open output file
276	if (fNameResult.IsNull() && fStoreResult)
277	{
278	TFile file = (TFile)gROOT->GetListOfFiles()->FindObject(fname);
279	if (file)
280	{
281	file->cd();
282	return WriteContainer(arr);
283	}
284	}
285
286	return WriteContainer(arr, fname, "UPDATE");
287	}
288
289	// --------------------------------------------------------------------------
290	//
291	// MJCut allows to setup several option by a resource file:
292	// MJCut.WriteSummary: yes, no
293	// MJCut.SummaryFile: filename
294	// MJCut.WriteResult: yes, no
295	// MJCut.ResultFile: filename
296	// MJCut.HistName: MHAlpha
297	//
298	Bool_t MJCut::CheckEnvLocal()
299	{
300	const TString f0(GetEnv("SummaryFile", ""));
301	const TString f1(GetEnv("ResultFile", ""));
302	if (!f0.IsNull())
303	SetNameSummaryFile(f0);
304	if (!f1.IsNull())
305	SetNameResultFile(f1);
306
307	EnableStorageOfSummary(GetEnv("SummaryFile", fStoreSummary));
308	EnableStorageOfResult(GetEnv("ResultFile", fStoreResult));
309	EnableWobbleMode(GetEnv("WobbleMode", fIsWobble));
310	EnableMonteCarloMode(GetEnv("MonteCarlo", fIsMonteCarlo));
311	EnableFullDisplay(GetEnv("FullDisplay", fFullDisplay));
312	//EnableSubstraction(GetEnv("HistogramSubstraction", fSubstraction));
313
314	SetNameHist(GetEnv("NameHist", fNameHist));
315	SetNameHistFS(GetEnv("NameHistFS", fNameHistFS));
316
317	return kTRUE;
318	}
319
320	// --------------------------------------------------------------------------
321	//
322	// Setup write to write:
323	// container tree optional?
324	// -------------- ---------- -----------
325	// "MHillas" to "Events"
326	// "MHillasSrc" to "Events"
327	// "MHadronness" to "Events" yes
328	// "MEnergyEst" to "Events" yes
329	// "DataType" to "Events"
330	//
331	void MJCut::SetupWriter(MWriteRootFile write, const char name) const
332	{
333	if (!write)
334	return;
335
336	write->SetName(name);
337	write->AddContainer("MHillas", "Events");
338	write->AddContainer("MHillasSrc", "Events");
339	write->AddContainer("MHillasExt", "Events");
340	write->AddContainer("MPointingPos", "Events");
341	write->AddContainer("MHillasSrcAnti", "Events", kFALSE);
342	write->AddContainer("MImagePar", "Events", kFALSE);
343	write->AddContainer("MNewImagePar", "Events", kFALSE);
344	write->AddContainer("MNewImagePar2", "Events", kFALSE);
345	write->AddContainer("MHadronness", "Events", kFALSE);
346	write->AddContainer("MSrcPosCam", "Events", kFALSE);
347	write->AddContainer("MSrcPosAnti", "Events", kFALSE);
348	write->AddContainer("ThetaSquared", "Events", kFALSE);
349	write->AddContainer("OpticalAxis", "Events", kFALSE);
350	write->AddContainer("Disp", "Events", kFALSE);
351	write->AddContainer("MTime", "Events", kFALSE);
352	write->AddContainer("MMcEvt", "Events", kFALSE);
353	write->AddContainer("DataType", "Events");
354	}
355
356	// --------------------------------------------------------------------------
357	//
358	// Create a new instance of an object with name name of class
359	// type fNameHist in parlist. It must derive from MHAlpha.
360	// Call ForceUsingSize for it and return its pointer.
361	// If something fails NULL is returned.
362	//
363	MHAlpha MJCut::CreateNewHist(MParList &plist, const char name) const
364	{
365	TClass *cls = gROOT->GetClass(fNameHist);
366	if (!cls)
367	{
368	*fLog << err << "Class " << fNameHist << " not found in dictionary... abort." << endl;
369	return NULL;
370	}
371	if (!cls->InheritsFrom(MHAlpha::Class()))
372	{
373	*fLog << err << "Class " << fNameHist << " doesn't inherit from MHAlpha... abort." << endl;
374	return NULL;
375	}
376
377	const TString objname(Form("Hist%s", name));
378	MHAlpha h = (MHAlpha)plist.FindCreateObj(fNameHist, objname);
379	if (!h)
380	return NULL;
381
382	h->ForceUsingSize();
383
384	return h;
385	}
386
387	// --------------------------------------------------------------------------
388	//
389	// Create a new instance of an object with name name of class
390	// type fNameHistFS in parlist. It must derive from MHFalseSource
391	// If something fails NULL is returned.
392	//
393	MH MJCut::CreateNewHistFS(MParList &plist, const char name) const
394	{
395	const TString cname(fNameHistFS.IsNull()?"MHFalseSource":fNameHistFS);
396
397	TClass *cls = gROOT->GetClass(cname);
398	if (!cls)
399	{
400	*fLog << err << "Class " << cname << " not found in dictionary... abort." << endl;
401	return NULL;
402	}
403	if (!cls->InheritsFrom("MHFalseSource"))
404	{
405	*fLog << err << "Class " << cname << " doesn't inherit from MHFalseSource... abort." << endl;
406	return NULL;
407	}
408
409	const TString objname(Form("FS%s", name));
410	return (MH*)plist.FindCreateObj(cname, objname);
411	}
412
413	Bool_t MJCut::Process(const MDataSet &set)
414	{
415	if (!set.IsValid())
416	{
417	*fLog << err << "ERROR - DataSet invalid!" << endl;
418	return kFALSE;
419	}
420
421	CheckEnv();
422
423	// --------------------------------------------------------------------------------
424
425	*fLog << inf;
426	fLog->Separator(GetDescriptor());
427	*fLog << "Perform cuts for data set " << set.GetName() << endl;
428	*fLog << endl;
429
430	// --------------------------------------------------------------------------------
431
432	// Setup Parlist
433	MParList plist;
434	plist.AddToList(this); // take care of fDisplay!
435
436	MParameterI par("DataType");
437	plist.AddToList(&par);
438
439	// Setup Tasklist
440	MTaskList tlist;
441	plist.AddToList(&tlist);
442
443	// La Palma Magic1
444	MObservatory obs;
445	plist.AddToList(&obs);
446
447	// Possible source position (eg. Wobble Mode)
448	MPointingPos source("MSourcePos");
449	if (set.HasSource())
450	{
451	if (!set.GetSourcePos(source))
452	return kFALSE;
453	plist.AddToList(&source);
454	*fLog << inf << "Using Source Position: ";
455	source.Print("RaDec");
456	}
457	else
458	*fLog << inf << "No source position applied..." << endl;
459
460	// Initialize default binnings
461	MBinning bins1(18, 0, 90, "BinningAlpha", "lin");
462	MBinning bins2(15, 10, 1e6 , "BinningSize", "log");
463	MBinning bins3(67, -0.005, 0.665, "BinningTheta", "asin");
464	MBinning bins4("BinningFalseSource");
465	MBinning bins5("BinningWidth");
466	MBinning bins6("BinningLength");
467	MBinning bins7("BinningDist");
468	MBinning bins8("BinningMaxDist");
469	MBinning bins9("BinningM3Long");
470	MBinning bins0("BinningConc1");
471	plist.AddToList(&bins1);
472	plist.AddToList(&bins2);
473	plist.AddToList(&bins3);
474	plist.AddToList(&bins4);
475	plist.AddToList(&bins5);
476	plist.AddToList(&bins6);
477	plist.AddToList(&bins7);
478	plist.AddToList(&bins8);
479	plist.AddToList(&bins9);
480	plist.AddToList(&bins0);
481	//plist.AddToList(&binsa);
482
483	// --------------------------------------------------------------------------------
484
485	// Setup fitter and histograms
486	MAlphaFitter fit;
487	plist.AddToList(&fit);
488	if (fIsWobble)
489	fit.SetScaleMode(MAlphaFitter::kNone);
490
491	MHAlpha *halphaoff = CreateNewHist(plist, "Off");
492	MFillH falpha(halphaoff, "", "FillHist");
493	MH *hfsoff = CreateNewHistFS(plist, "Off");
494	MFillH ffs(hfsoff, "MHillas", "FillFS");
495
496	// FIXME: If fPathIn read cuts and energy estimator from file!
497	MContinue cont0("", "Cut0");
498	MContinue cont1("", "Cut1");
499	MContinue cont2("", "Cut2");
500	MContinue cont3("", "Cut3");
501	cont0.SetAllowEmpty();
502	cont1.SetAllowEmpty();
503	cont2.SetAllowEmpty();
504	cont3.SetAllowEmpty();
505
506	// ------------- Loop Off Data --------------------
507	MReadReports readoffdata;
508	readoffdata.AddTree("Events", "MTime.", MReadReports::kMaster);
509	readoffdata.AddTree("Drive", MReadReports::kRequired);
510	readoffdata.AddTree("Starguider", MReadReports::kRequired);
511	readoffdata.AddTree("EffectiveOnTime");
512
513	MReadMarsFile readoffmc("Events");
514	readoffmc.DisableAutoScheme();
515
516	MRead &readoff = fIsMonteCarlo ? (MRead&)readoffmc : (MRead&)readoffdata;
517	if (fIsWobble)
518	set.AddFilesOn(readoff);
519	else
520	set.AddFilesOff(readoff);
521
522	const TString path(Form("%s/", fPathOut.Data()));
523	TString fname0(path);
524	TString fname1(path);
525	fname0 += fNameSummary.IsNull() ? (TString) Form("ganymed%08d-summary.root", set.GetNumAnalysis()) : fNameSummary;
526	fname1 += fNameResult.IsNull() ? (TString) Form("ganymed%08d.root", set.GetNumAnalysis()) : fNameResult;
527
528	MWriteRootFile *write0 = CanStoreSummary() ? new MWriteRootFile(fPathOut.IsNull()?0:fname0.Data(), fOverwrite?"RECREATE":"NEW") : 0;
529	MWriteRootFile *write1 = CanStoreResult() ? new MWriteRootFile(fPathOut.IsNull()?0:fname1.Data(), fOverwrite?"RECREATE":"NEW") : 0;
530	SetupWriter(write0, "WriteAfterCut0");
531	SetupWriter(write1, "WriteAfterCut3");
532
533	/*
534	MEnergyEstimate est;
535
536	MTaskEnv taskenv1("EstimateEnergy");
537	taskenv1.SetDefault(fEstimateEnergy ? fEstimateEnergy : &est);
538	*/
539	MTaskEnv taskenv2("CalcHadronness");
540	taskenv2.SetDefault(fCalcHadronness);
541
542	MFillH fill1a("MHHillasOffPre [MHHillas]", "MHillas", "FillHillasPre");
543	MFillH fill2a("MHHillasOffPost [MHHillas]", "MHillas", "FillHillasPost");
544	MFillH fill3a("MHVsSizeOffPost [MHVsSize]", "MHillasSrc", "FillVsSizePost");
545	MFillH fill4a("MHHilExtOffPost [MHHillasExt]", "MHillasSrc", "FillHilExtPost");
546	MFillH fill5a("MHHilSrcOffPost [MHHillasSrc]", "MHillasSrc", "FillHilSrcPost");
547	MFillH fill6a("MHImgParOffPost [MHImagePar]", "MImagePar", "FillImgParPost");
548	MFillH fill7a("MHNewParOffPost [MHNewImagePar]", "MNewImagePar", "FillNewParPost");
549	fill1a.SetNameTab("PreCut");
550	fill2a.SetNameTab("PostCut");
551	fill3a.SetNameTab("VsSize");
552	fill4a.SetNameTab("HilExt");
553	fill5a.SetNameTab("HilSrc");
554	fill6a.SetNameTab("ImgPar");
555	fill7a.SetNameTab("NewPar");
556
557	MPrint print2("MEffectiveOnTime");
558	print2.EnableSkip();
559
560	// How to get source position from off- and on-data?
561	MSrcPosCalc scalc;
562	scalc.SetMode(fIsWobble?MSrcPosCalc::kWobble:MSrcPosCalc::kOffData); /********************/
563
564	MSrcPosCorrect scor;
565
566	MHillasCalc hcalc;
567	MHillasCalc hcalc2("MHillasCalcAnti");
568	hcalc.SetFlags(MHillasCalc::kCalcHillasSrc);
569	hcalc2.SetFlags(MHillasCalc::kCalcHillasSrc);
570	hcalc2.SetNameHillasSrc("MHillasSrcAnti");
571	hcalc2.SetNameSrcPosCam("MSrcPosAnti");
572
573	MTaskList tlist2;
574	tlist2.AddToList(&scalc);
575	tlist2.AddToList(&scor);
576	tlist2.AddToList(&hcalc);
577	if (fIsWobble)
578	tlist2.AddToList(&hcalc2);
579	//tlist2.AddToList(&taskenv1);
580	tlist2.AddToList(&taskenv2);
581	tlist2.AddToList(&cont0);
582	if (write0)
583	tlist2.AddToList(write0);
584	if (!fWriteOnly)
585	tlist2.AddToList(&fill1a);
586	tlist2.AddToList(&cont1);
587	if (!fWriteOnly && (!fIsWobble \|\| !fNameHistFS.IsNull()))
588	tlist2.AddToList(&ffs);
589	tlist2.AddToList(&cont2);
590	if (!fWriteOnly)
591	{
592	tlist2.AddToList(&fill2a);
593	if (fFullDisplay)
594	{
595	tlist2.AddToList(&fill3a);
596	tlist2.AddToList(&fill4a);
597	tlist2.AddToList(&fill5a);
598	tlist2.AddToList(&fill6a);
599	tlist2.AddToList(&fill7a);
600	}
601	}
602	if (!fWriteOnly)
603	tlist2.AddToList(&falpha);
604	tlist2.AddToList(&cont3);
605	if (write1)
606	tlist2.AddToList(write1);
607
608	MPointingDevCalc devcalc;
609
610	tlist.AddToList(&readoff);
611	if (gLog.GetDebugLevel()>4)
612	tlist.AddToList(&print2, "EffectiveOnTime");
613	tlist.AddToList(&devcalc, "Starguider");
614	tlist.AddToList(&tlist2, "Events");
615
616	par.SetVal(0);
617
618	// Create and setup the eventloop
619	MEvtLoop evtloop(fName);
620	evtloop.SetParList(&plist);
621	evtloop.SetDisplay(fDisplay);
622	evtloop.SetLogStream(fLog);
623	if (!SetupEnv(evtloop))
624	return kFALSE;
625
626	TObjArray cont;
627	cont.Add(&cont0);
628	cont.Add(&cont1);
629	cont.Add(&cont2);
630	cont.Add(&cont3);
631	//if (taskenv1.GetTask())
632	// cont.Add(taskenv1.GetTask());
633	if (taskenv2.GetTask())
634	cont.Add(taskenv2.GetTask());
635
636	if (!WriteTasks(set.GetNumAnalysis(), cont))
637	return kFALSE;
638
639	if (set.HasOffSequences() \|\| fIsWobble)
640	{
641	// Execute first analysis
642	if (!evtloop.Eventloop(fMaxEvents))
643	{
644	*fLog << err << GetDescriptor() << ": Processing of off-sequences failed." << endl;
645	return kFALSE;
646	}
647
648	if (!evtloop.GetDisplay())
649	{
650	*fLog << err << GetDescriptor() << ": Execution stopped by user." << endl;
651	return kFALSE;
652	}
653
654	//plist.FindObject("MTimeEffectiveOnTime")->Clear();
655	}
656	else
657	{
658	// This is the simplest way to remove the two object from the parlist
659	delete halphaoff;
660	delete hfsoff;
661	}
662
663	// ------------- Loop On Data --------------------
664	MReadReports readondata;
665	readondata.AddTree("Events", "MTime.", MReadReports::kMaster);
666	readondata.AddTree("Drive", MReadReports::kRequired);
667	readondata.AddTree("Starguider", MReadReports::kRequired);
668	readondata.AddTree("EffectiveOnTime");
669
670	MReadMarsFile readonmc("Events");
671	readonmc.DisableAutoScheme();
672
673	MRead &readon = fIsMonteCarlo ? (MRead&)readonmc : (MRead&)readondata;
674	set.AddFilesOn(readon);
675
676	scalc.SetMode(MSrcPosCalc::kDefault);
677
678	MFillH fill1b("MHHillasOnPre [MHHillas]", "MHillas", "FillHillasPre");
679	MFillH fill2b("MHHillasOnPost [MHHillas]", "MHillas", "FillHillasPost");
680	MFillH fill3b("MHVsSizeOnPost [MHVsSize]", "MHillasSrc", "FillVsSizePost");
681	MFillH fill4b("MHHilExtOnPost [MHHillasExt]", "MHillasSrc", "FillHilExtPost");
682	MFillH fill5b("MHHilSrcOnPost [MHHillasSrc]", "MHillasSrc", "FillHilSrcPost");
683	MFillH fill6b("MHImgParOnPost [MHImagePar]", "MImagePar", "FillImgParPost");
684	MFillH fill7b("MHNewParOnPost [MHNewImagePar]", "MNewImagePar", "FillNewParPost");
685	fill1b.SetNameTab("PreCut");
686	fill2b.SetNameTab("PostCut");
687	fill3b.SetNameTab("VsSize");
688	fill4b.SetNameTab("HilExt");
689	fill5b.SetNameTab("HilSrc");
690	fill6b.SetNameTab("ImgPar");
691	fill7b.SetNameTab("NewPar");
692	fill1b.SetDrawOption(set.HasOffSequences()\|\|fIsWobble?"same":"");
693	fill2b.SetDrawOption(set.HasOffSequences()\|\|fIsWobble?"same":"");
694	fill3b.SetDrawOption(set.HasOffSequences()\|\|fIsWobble?"same":"");
695	fill4b.SetDrawOption(set.HasOffSequences()\|\|fIsWobble?"same":"");
696	fill5b.SetDrawOption(set.HasOffSequences()\|\|fIsWobble?"same":"");
697	fill6b.SetDrawOption(set.HasOffSequences()\|\|fIsWobble?"same":"");
698	fill7b.SetDrawOption(set.HasOffSequences()\|\|fIsWobble?"same":"");
699
700	/*
701	MHVsTime hvs("MEffectiveOnTime.fVal");
702	hvs.SetTitle("Effective On-Time vs. Time;;T_{on}");
703	MFillH fillvs(&hvs, "MTimeEffectiveOnTime", "FillOnTime");
704	fillvs.SetNameTab("OnTime");
705	*/
706	MH3 hvs("MPointingPos.fZd");
707	hvs.SetName("Theta");
708	hvs.SetTitle("Effective On-Time vs. Zenith Angle;\\Theta [\\circ];T_{on} [s]");
709	MFillH fillvs(&hvs, "", "FillOnTime");
710	if (!fIsMonteCarlo)
711	fillvs.SetWeight("MEffectiveOnTime");
712	fillvs.SetNameTab("OnTime");
713
714	/*
715	MParameterD weight;
716	weight.SetVal(-1);
717	fill2a.SetWeight(&weight);
718	fill3a.SetWeight(&weight);
719	fill4a.SetWeight(&weight);
720	fill5a.SetWeight(&weight);
721	fill6a.SetWeight(&weight);
722	fill7a.SetWeight(&weight);
723	if (fSubstraction)
724	{
725	fill2a.SetNameTab("PostCut-");
726	fill3a.SetNameTab("VsSize-");
727	fill4a.SetNameTab("HilExt-");
728	fill5a.SetNameTab("HilSrc-");
729	fill6a.SetNameTab("ImgPar-");
730	fill7a.SetNameTab("NewPar-");
731	}
732	*/
733	MHAlpha *halphaon=CreateNewHist(plist);
734	MFillH falpha2(halphaon, "", "FillHist");
735	MH *hfs=CreateNewHistFS(plist);
736	MFillH ffs2(hfs, "MHillas", "FillFS");
737	MFillH fillphi("MHPhi", "", "FillPhi");
738	fillphi.SetDrawOption("anticut");
739
740	tlist.Replace(&readon);
741	if (!fWriteOnly)
742	{
743	tlist2.Replace(&fill1b);
744	/* if (fIsWobble)
745	{
746	tlist2.AddToListAfter(&fill2b, &fill2a);
747	tlist2.AddToListAfter(&fill3b, &fill3a);
748	}
749	else
750	*/
751	tlist2.Replace(&fill2b);
752	if (fFullDisplay)
753	{
754	tlist2.Replace(&fill3b);
755	tlist2.Replace(&fill4b);
756	tlist2.Replace(&fill5b);
757	tlist2.Replace(&fill6b);
758	tlist2.Replace(&fill7b);
759	}
760	tlist2.Replace(&falpha2);
761	if (!fIsWobble \|\| !fNameHist.IsNull())
762	tlist2.Replace(&ffs2);
763	if (fIsWobble)
764	{
765	tlist2.AddToListAfter(&fillphi, &falpha2);
766	if (!fNameHist.IsNull())
767	tlist2.RemoveFromList(&ffs);
768	}
769
770	if (!fIsMonteCarlo)
771	tlist.AddToList(&fillvs, "EffectiveOnTime");
772	else
773	tlist2.AddToListBefore(&fillvs, &scalc);
774	}
775
776	par.SetVal(1);
777
778	// Execute first analysis
779	if (!evtloop.Eventloop(fMaxEvents))
780	{
781	*fLog << err << GetDescriptor() << ": Processing of on-sequences failed." << endl;
782	return kFALSE;
783	}
784
785	if (write0)
786	delete write0;
787	if (write1)
788	delete write1;
789
790	// FIXME: Perform fit and plot energy dependant alpha plots
791	// and fit result to new tabs!
792	if (!WriteResult(plist, set.GetNumAnalysis()))
793	return kFALSE;
794
795	*fLog << all << GetDescriptor() << ": Done." << endl;
796	*fLog << endl << endl;
797
798	return kTRUE;
799	}

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

Download in other formats: