Context Navigation

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

Visit:

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

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

Download in other formats: