Context Navigation

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

Visit:

Last change on this file since 8244 was 8244, checked in by tbretz, 18 years ago
* empty log message *
File size: 29.7 KB

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

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

Download in other formats: