Context Navigation

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

Visit:

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

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

Download in other formats: