Context Navigation

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

Visit:

Last change on this file since 8441 was 8441, checked in by tbretz, 17 years ago
* empty log message *
File size: 30.4 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), 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	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("RunNumber", "Events");
378	write->AddContainer("EvtNumber", "Events");
379	// write->AddContainer("MMuonSearchPar", "Events", kFALSE);
380	// write->AddContainer("MMuonCalibPar", "Events", kFALSE);
381	}
382
383	// --------------------------------------------------------------------------
384	//
385	// Create a new instance of an object with name name of class
386	// type fNameHist in parlist. It must derive from MHAlpha.
387	// Call ForceUsingSize for it and return its pointer.
388	// If something fails NULL is returned.
389	//
390	MHAlpha MJCut::CreateNewHist(MParList &plist, const char name) const
391	{
392	TClass *cls = gROOT->GetClass(fNameHist);
393	if (!cls)
394	{
395	*fLog << err << "Class " << fNameHist << " not found in dictionary... abort." << endl;
396	return NULL;
397	}
398	if (!cls->InheritsFrom(MHAlpha::Class()))
399	{
400	*fLog << err << "Class " << fNameHist << " doesn't inherit from MHAlpha... abort." << endl;
401	return NULL;
402	}
403
404	const TString objname(Form("Hist%s", name));
405	MHAlpha h = (MHAlpha)plist.FindCreateObj(fNameHist, objname);
406	if (!h)
407	return NULL;
408
409	h->ForceUsingSize();
410
411	return h;
412	}
413
414	// --------------------------------------------------------------------------
415	//
416	// Create a new instance of an object with name name of class
417	// type fNameHistFS in parlist. It must derive from MHFalseSource
418	// If something fails NULL is returned.
419	//
420	MH MJCut::CreateNewHistFS(MParList &plist, const char name) const
421	{
422	const TString cname(fNameHistFS.IsNull()?"MHFalseSource":fNameHistFS.Data());
423
424	TClass *cls = gROOT->GetClass(cname);
425	if (!cls)
426	{
427	*fLog << err << "Class " << cname << " not found in dictionary... abort." << endl;
428	return NULL;
429	}
430	if (!cls->InheritsFrom("MHFalseSource"))
431	{
432	*fLog << err << "Class " << cname << " doesn't inherit from MHFalseSource... abort." << endl;
433	return NULL;
434	}
435
436	const TString objname(Form("FS%s", name));
437	return (MH*)plist.FindCreateObj(cname, objname);
438	}
439
440	Bool_t MJCut::FillSrcPosCam(const MDataSet &set, MPointingPos &source, MHSrcPosCam &hsrcpos)
441	{
442	*fLog << inf;
443	fLog->Separator(GetDescriptor());
444	*fLog << "Filling MHSrcPosCam " << set.GetName() << endl;
445	*fLog << endl;
446
447	// --------------------------------------------------------------------------------
448
449	// Setup Parlist
450	MParList plist;
451	plist.AddToList(this); // take care of fDisplay!
452
453	// Setup Tasklist
454	MTaskList tlist;
455	plist.AddToList(&tlist);
456
457	// La Palma Magic1, Possible source position
458	MObservatory obs;
459	plist.AddToList(&obs);
460	plist.AddToList(&source);
461
462	// Initialize default binnings
463	// MBinning bins1(18, 0, 90, "BinningSrcPosCam", "lin");
464	// plist.AddToList(&bins1);
465
466	// ------------- Loop Off Data --------------------
467	MReadReports read;
468
469	read.EnableAutoScheme();
470	read.AddToBranchList("MTimeEffectiveOnTime.*");
471	read.AddToBranchList("MEffectiveOnTime.*");
472
473	read.AddTree("Events", "MTime.", MReadReports::kMaster);
474	read.AddTree("Drive", MReadReports::kRequired);
475	read.AddTree("Starguider", MReadReports::kRequired);
476	read.AddTree("EffectiveOnTime");
477
478	if (!set.AddFilesOn(read))
479	return kFALSE;
480
481	MFillH fill(&hsrcpos, "MSrcPosCam", "FillSrcPosCam");
482	fill.SetNameTab("SrcPos");
483
484	// How to get source position from off- and on-data?
485	MSrcPosCorrect scor;
486	MSrcPosCalc scalc;
487	scalc.SetMode(MSrcPosCalc::kDefault);
488
489	MPointingDevCalc devcalc;
490
491	tlist.AddToList(&read);
492	tlist.AddToList(&devcalc, "Starguider");
493	tlist.AddToList(&scalc, "Events");
494	tlist.AddToList(&scor, "Events");
495	tlist.AddToList(&fill, "Events");
496
497	// by setting it here it is distributed to all consecutive tasks
498	tlist.SetAccelerator(MTask::kAccDontReset\|MTask::kAccDontTime);
499
500	// Create and setup the eventloop
501	MEvtLoop evtloop(fName);
502	evtloop.SetParList(&plist);
503	evtloop.SetDisplay(fDisplay);
504	evtloop.SetLogStream(fLog);
505	if (!SetupEnv(evtloop))
506	return kFALSE;
507
508	// Execute first analysis
509	if (!evtloop.Eventloop(fMaxEvents))
510	{
511	*fLog << err << GetDescriptor() << ": Processing of on-sequences to fill SrcPosCam failed." << endl;
512	return kFALSE;
513	}
514
515	if (!evtloop.GetDisplay())
516	{
517	*fLog << err << GetDescriptor() << ": Execution stopped by user." << endl;
518	return kFALSE;
519	}
520
521	*fLog << all << GetDescriptor() << ": Done." << endl;
522	*fLog << endl << endl;
523
524	return kTRUE;
525	}
526
527	Int_t MJCut::Process(const MDataSet &set)
528	{
529	if (!set.IsValid())
530	{
531	*fLog << err << "ERROR - DataSet invalid!" << endl;
532	return kFALSE;
533	}
534
535	CheckEnv();
536
537	// --------------------------------------------------------------------------------
538
539	// Possible source position (eg. Wobble Mode)
540	MPointingPos source("MSourcePos");
541	if (set.HasSource())
542	{
543	if (!set.GetSourcePos(source))
544	return -1;
545	*fLog << all;
546	source.Print("RaDec");
547	}
548	else
549	*fLog << all << "No source position applied..." << endl;
550
551	MParList plist;
552
553	MHSrcPosCam hsrcpos(fIsWobble);
554	if (!fIsWobble && source.IsInitialized() && fRndmSrcPos)
555	{
556	if (!FillSrcPosCam(set, source, hsrcpos))
557	return -2;
558	plist.AddToList(&hsrcpos);
559	}
560
561	// --------------------------------------------------------------------------------
562
563	*fLog << inf;
564	fLog->Separator(GetDescriptor());
565	*fLog << "Perform cuts for data set " << set.GetName() << endl;
566	*fLog << endl;
567
568	// --------------------------------------------------------------------------------
569
570	// Setup Parlist
571	plist.AddToList(this); // take care of fDisplay!
572
573	MParameterI par("DataType");
574	plist.AddToList(&par);
575
576	// Setup Tasklist
577	MTaskList tlist;
578	plist.AddToList(&tlist);
579
580	// La Palma Magic1
581	MObservatory obs;
582	plist.AddToList(&obs);
583
584	if (source.IsInitialized())
585	plist.AddToList(&source);
586
587	// Initialize default binnings
588	MBinning bins1( 18, 0, 90, "BinningAlpha", "lin");
589	MBinning bins2( 15, 10, 1e6 , "BinningSize", "log");
590	MBinning bins3( 67, -0.005, 0.665, "BinningTheta", "asin");
591	//MBinning binsT(150, 0, 150, "BinningDeltaT", "lin");
592	MBinning bins4("BinningFalseSource");
593	MBinning bins5("BinningWidth");
594	MBinning bins6("BinningLength");
595	MBinning bins7("BinningDist");
596	MBinning bins8("BinningMaxDist");
597	MBinning bins9("BinningM3Long");
598	MBinning bins0("BinningConc1");
599	plist.AddToList(&bins1);
600	plist.AddToList(&bins2);
601	plist.AddToList(&bins3);
602	plist.AddToList(&bins4);
603	plist.AddToList(&bins5);
604	plist.AddToList(&bins6);
605	plist.AddToList(&bins7);
606	plist.AddToList(&bins8);
607	plist.AddToList(&bins9);
608	plist.AddToList(&bins0);
609	//plist.AddToList(&binsT);
610
611	// --------------------------------------------------------------------------------
612
613	// Setup fitter and histograms
614	MAlphaFitter fit;
615	plist.AddToList(&fit);
616	if (fIsWobble)
617	fit.SetScaleMode(MAlphaFitter::kNone);
618
619	MHAlpha *halphaoff = CreateNewHist(plist, "Off");
620	MFillH falpha(halphaoff, "", "FillHist");
621	MH *hfsoff = CreateNewHistFS(plist, "Off");
622	MFillH ffs(hfsoff, "MHillas", "FillFS");
623
624	// FIXME: If fPathIn read cuts and energy estimator from file!
625	MContinue cont0("", "Cut0");
626	MContinue cont1("", "Cut1");
627	MContinue cont2("", "Cut2");
628	MContinue cont3("", "Cut3");
629	cont0.SetAllowEmpty();
630	cont1.SetAllowEmpty();
631	cont2.SetAllowEmpty();
632	cont3.SetAllowEmpty();
633
634	// ------------- Loop Off Data --------------------
635	MReadReports readoffdata;
636	readoffdata.AddTree("Events", "MTime.", MReadReports::kMaster);
637	readoffdata.AddTree("Drive", MReadReports::kRequired);
638	readoffdata.AddTree("Starguider", MReadReports::kRequired);
639	readoffdata.AddTree("EffectiveOnTime");
640
641	MReadMarsFile readoffmc("Events");
642	readoffmc.DisableAutoScheme();
643
644	MRead &readoff = set.IsMonteCarlo() ? (MRead&)readoffmc : (MRead&)readoffdata;
645	const Bool_t setrc = fIsWobble ? set.AddFilesOn(readoff) : set.AddFilesOff(readoff);
646	if (!setrc && set.HasOffSequences())
647	{
648	*fLog << err << "MDataSet::AddFiles" << (fIsWobble?"On":"Off") << " failed." << endl;
649	return kFALSE;
650	}
651
652	const TString path(Form("%s/", fPathOut.Data()));
653	TString fname0(path);
654	TString fname1(path);
655	fname0 += fNameSummary.IsNull() ? (TString) Form("ganymed%08d-summary.root", set.GetNumAnalysis()) : fNameSummary;
656	fname1 += fNameResult.IsNull() ? (TString) Form("ganymed%08d.root", set.GetNumAnalysis()) : fNameResult;
657
658	MWriteRootFile *write0 = CanStoreSummary() ? new MWriteRootFile(fPathOut.IsNull()?0:fname0.Data(), fOverwrite?"RECREATE":"NEW") : 0;
659	MWriteRootFile *write1 = CanStoreResult() ? new MWriteRootFile(fPathOut.IsNull()?0:fname1.Data(), fOverwrite?"RECREATE":"NEW") : 0;
660	SetupWriter(write0, "WriteAfterCut0");
661	SetupWriter(write1, "WriteAfterCut3");
662
663	MTaskEnv taskenv2("CalcHadronness");
664	taskenv2.SetDefault(fCalcHadronness);
665
666	MTaskEnv taskenv3("CalcDisp");
667	taskenv3.SetDefault(fCalcDisp);
668
669	MParameterCalc setevtnum("MRawEvtHeader.fDAQEvtNumber", "SetEvtNumber");
670	setevtnum.SetNameParameter("EvtNumber");
671
672	MParameterCalc setrunnum("MRawRunHeader.fRunNumber", "SetRunNumber");
673	setrunnum.SetNameParameter("RunNumber");
674
675	// MFillH fill0a("OnPos [MHSrcPosCam]", "MSrcPosCam", "FillSrcPosCam");
676	MFillH fill1a("MHHillasOffPre [MHHillas]", "MHillas", "FillHillasPre");
677	MFillH fill2a("MHHillasOffPost [MHHillas]", "MHillas", "FillHillasPost");
678	MFillH fill3a("MHVsSizeOffPost [MHVsSize]", "MHillasSrc", "FillVsSizePost");
679	MFillH fill4a("MHHilExtOffPost [MHHillasExt]", "MHillasSrc", "FillHilExtPost");
680	MFillH fill5a("MHHilSrcOffPost [MHHillasSrc]", "MHillasSrc", "FillHilSrcPost");
681	MFillH fill6a("MHImgParOffPost [MHImagePar]", "MImagePar", "FillImgParPost");
682	MFillH fill7a("MHNewParOffPost [MHNewImagePar]", "MNewImagePar", "FillNewParPost");
683	//MFillH fill9a("MHEffOffTime [MHEffectiveOnTime]", "MTime", "FillEffOnTime");
684	fill1a.SetNameTab("PreCut");
685	fill2a.SetNameTab("PostCut");
686	fill3a.SetNameTab("VsSize");
687	fill4a.SetNameTab("HilExt");
688	fill5a.SetNameTab("HilSrc");
689	fill6a.SetNameTab("ImgPar");
690	fill7a.SetNameTab("NewPar");
691	//fill9a.SetNameTab("EffOffT");
692
693	//MFDataMember fbin("Bin.fVal", '>', 0);
694	//fill9a.SetFilter(&fbin);
695
696	MPrint print2("MEffectiveOnTime");
697	print2.EnableSkip();
698
699	// How to get source position from off- and on-data?
700	MSrcPosCalc scalc;
701	scalc.SetMode(fIsWobble?MSrcPosCalc::kWobble:MSrcPosCalc::kOffData); /********************/
702
703	MSrcPosCorrect scor;
704
705	MHillasCalc hcalc;
706	MHillasCalc hcalc2("MHillasCalcAnti");
707	hcalc.SetFlags(MHillasCalc::kCalcHillasSrc);
708	hcalc2.SetFlags(MHillasCalc::kCalcHillasSrc);
709	hcalc2.SetNameHillasSrc("MHillasSrcAnti");
710	hcalc2.SetNameSrcPosCam("MSrcPosAnti");
711
712	MSrcPosRndm srcrndm;
713
714	MTaskList tlist2;
715	tlist2.AddToList(&scalc);
716	tlist2.AddToList(&scor);
717	tlist2.AddToList(&srcrndm);
718	//if (fRndmSrcPos && !fIsWobble)
719	// tlist2.AddToList(&fill0a);
720	tlist2.AddToList(&hcalc);
721	if (fIsWobble)
722	tlist2.AddToList(&hcalc2);
723	//tlist2.AddToList(&taskenv1);
724	tlist2.AddToList(&cont0);
725	tlist2.AddToList(&taskenv2);
726	tlist2.AddToList(&taskenv3);
727	tlist2.AddToList(&setrunnum);
728	tlist2.AddToList(&setevtnum);
729	if (write0)
730	tlist2.AddToList(write0);
731	if (!fWriteOnly)
732	tlist2.AddToList(&fill1a);
733	tlist2.AddToList(&cont1);
734	if (!fWriteOnly && (!fIsWobble \|\| !fNameHistFS.IsNull()))
735	tlist2.AddToList(&ffs);
736	tlist2.AddToList(&cont2);
737	if (!fWriteOnly)
738	{
739	tlist2.AddToList(&fill2a);
740	if (fFullDisplay)
741	{
742	tlist2.AddToList(&fill3a);
743	tlist2.AddToList(&fill4a);
744	tlist2.AddToList(&fill5a);
745	tlist2.AddToList(&fill6a);
746	tlist2.AddToList(&fill7a);
747	}
748	}
749	if (!fWriteOnly)
750	{
751	tlist2.AddToList(&falpha);
752	/* if (!fIsMonteCarlo)
753	{
754	tlist2.AddToList(&fbin);
755	tlist2.AddToList(&fill9a);
756	} */
757	}
758	tlist2.AddToList(&cont3);
759	if (write1)
760	tlist2.AddToList(write1);
761
762	MPointingDevCalc devcalc;
763
764	tlist.AddToList(&readoff);
765	if (gLog.GetDebugLevel()>4)
766	tlist.AddToList(&print2, "EffectiveOnTime");
767	tlist.AddToList(&devcalc, "Starguider");
768	tlist.AddToList(&tlist2, "Events");
769
770	// by setting it here it is distributed to all consecutive tasks
771	tlist.SetAccelerator(MTask::kAccDontReset\|MTask::kAccDontTime);
772
773	par.SetVal(0);
774
775	// Create and setup the eventloop
776	MEvtLoop evtloop(fName);
777	evtloop.SetParList(&plist);
778	evtloop.SetDisplay(fDisplay);
779	evtloop.SetLogStream(fLog);
780	if (!SetupEnv(evtloop))
781	return kFALSE;
782
783	TObjArray cont;
784	cont.Add(&cont0);
785	cont.Add(&cont1);
786	cont.Add(&cont2);
787	cont.Add(&cont3);
788	//if (taskenv1.GetTask())
789	// cont.Add(taskenv1.GetTask());
790	if (taskenv2.GetTask())
791	cont.Add(taskenv2.GetTask());
792	if (taskenv3.GetTask())
793	cont.Add(taskenv3.GetTask());
794
795	if (!WriteTasks(set.GetNumAnalysis(), cont))
796	return kFALSE;
797
798	if (set.HasOffSequences() \|\| fIsWobble)
799	{
800	// Execute first analysis
801	if (!evtloop.Eventloop(fMaxEvents))
802	{
803	*fLog << err << GetDescriptor() << ": Processing of off-sequences failed." << endl;
804	return -3;
805	}
806
807	if (!evtloop.GetDisplay())
808	{
809	*fLog << err << GetDescriptor() << ": Execution stopped by user." << endl;
810	return kFALSE;
811	}
812
813	//plist.FindObject("MTimeEffectiveOnTime")->Clear();
814	}
815	else
816	{
817	// This is the simplest way to remove the two object from the parlist
818	delete halphaoff;
819	delete hfsoff;
820	}
821
822	// ------------- Loop On Data --------------------
823	MReadReports readondata;
824	readondata.AddTree("Events", "MTime.", MReadReports::kMaster);
825	readondata.AddTree("Drive", MReadReports::kRequired);
826	readondata.AddTree("Starguider", MReadReports::kRequired);
827	readondata.AddTree("EffectiveOnTime");
828
829	MReadMarsFile readonmc("Events");
830	readonmc.DisableAutoScheme();
831
832	MRead &readon = set.IsMonteCarlo() ? (MRead&)readonmc : (MRead&)readondata;
833	if (!set.AddFilesOn(readon))
834	return kFALSE;
835
836	scalc.SetMode(MSrcPosCalc::kDefault);
837
838	MFillH fill1b("MHHillasOnPre [MHHillas]", "MHillas", "FillHillasPre");
839	MFillH fill2b("MHHillasOnPost [MHHillas]", "MHillas", "FillHillasPost");
840	MFillH fill3b("MHVsSizeOnPost [MHVsSize]", "MHillasSrc", "FillVsSizePost");
841	MFillH fill4b("MHHilExtOnPost [MHHillasExt]", "MHillasSrc", "FillHilExtPost");
842	MFillH fill5b("MHHilSrcOnPost [MHHillasSrc]", "MHillasSrc", "FillHilSrcPost");
843	MFillH fill6b("MHImgParOnPost [MHImagePar]", "MImagePar", "FillImgParPost");
844	MFillH fill7b("MHNewParOnPost [MHNewImagePar]", "MNewImagePar", "FillNewParPost");
845	//MFillH fill9b("MHEffOnTime [MHEffectiveOnTime]", "MTime", "FillEffOnTime");
846	fill1b.SetNameTab("PreCut");
847	fill2b.SetNameTab("PostCut");
848	fill3b.SetNameTab("VsSize");
849	fill4b.SetNameTab("HilExt");
850	fill5b.SetNameTab("HilSrc");
851	fill6b.SetNameTab("ImgPar");
852	fill7b.SetNameTab("NewPar");
853	//fill9b.SetNameTab("EffOnT");
854	fill1b.SetDrawOption(set.HasOffSequences()\|\|fIsWobble?"same":"");
855	fill2b.SetDrawOption(set.HasOffSequences()\|\|fIsWobble?"same":"");
856	fill3b.SetDrawOption(set.HasOffSequences()\|\|fIsWobble?"same":"");
857	fill4b.SetDrawOption(set.HasOffSequences()\|\|fIsWobble?"same":"");
858	fill5b.SetDrawOption(set.HasOffSequences()\|\|fIsWobble?"same":"");
859	fill6b.SetDrawOption(set.HasOffSequences()\|\|fIsWobble?"same":"");
860	fill7b.SetDrawOption(set.HasOffSequences()\|\|fIsWobble?"same":"");
861
862	//fill9b.SetFilter(&fbin);
863
864	/*
865	MHVsTime hvs("MEffectiveOnTime.fVal");
866	hvs.SetTitle("Effective On-Time vs. Time;;T_{on}");
867	MFillH fillvs(&hvs, "MTimeEffectiveOnTime", "FillOnTime");
868	fillvs.SetNameTab("OnTime");
869	*/
870	MH3 hvs("MPointingPos.fZd");
871	hvs.SetName("Theta");
872	hvs.SetTitle("Effective On-Time vs. Zenith Angle;\\Theta [\\circ];T_{on} [s]");
873	MFillH fillvs(&hvs, "", "FillOnTime");
874	if (!set.IsMonteCarlo())
875	fillvs.SetWeight("MEffectiveOnTime");
876	fillvs.SetNameTab("OnTime");
877
878	/*
879	MParameterD weight;
880	weight.SetVal(-1);
881	fill2a.SetWeight(&weight);
882	fill3a.SetWeight(&weight);
883	fill4a.SetWeight(&weight);
884	fill5a.SetWeight(&weight);
885	fill6a.SetWeight(&weight);
886	fill7a.SetWeight(&weight);
887	if (fSubstraction)
888	{
889	fill2a.SetNameTab("PostCut-");
890	fill3a.SetNameTab("VsSize-");
891	fill4a.SetNameTab("HilExt-");
892	fill5a.SetNameTab("HilSrc-");
893	fill6a.SetNameTab("ImgPar-");
894	fill7a.SetNameTab("NewPar-");
895	}
896	*/
897	MHAlpha *halphaon=CreateNewHist(plist);
898	MFillH falpha2(halphaon, "", "FillHist");
899	MH *hfs=CreateNewHistFS(plist);
900	MFillH ffs2(hfs, "MHillas", "FillFS");
901	MFillH fillphi("MHPhi", "", "FillPhi");
902	fillphi.SetDrawOption("anticut");
903
904	tlist.Replace(&readon);
905	if (fRndmSrcPos && !fIsWobble)
906	{
907	// tlist2.RemoveFromList(&fill0a);
908	tlist2.RemoveFromList(&srcrndm);
909	}
910
911	MFillH fillsrc(&hsrcpos, "MSrcPosCam", "FillSrcPosCam");
912	fillsrc.SetNameTab("SrcPos");
913	if (fIsWobble)
914	tlist2.AddToListBefore(&fillsrc, &hcalc);
915
916	if (!fWriteOnly)
917	{
918	tlist2.Replace(&fill1b);
919	/* if (fIsWobble)
920	{
921	tlist2.AddToListAfter(&fill2b, &fill2a);
922	tlist2.AddToListAfter(&fill3b, &fill3a);
923	}
924	else
925	*/
926	tlist2.Replace(&fill2b);
927	if (fFullDisplay)
928	{
929	tlist2.Replace(&fill3b);
930	tlist2.Replace(&fill4b);
931	tlist2.Replace(&fill5b);
932	tlist2.Replace(&fill6b);
933	tlist2.Replace(&fill7b);
934	}
935	tlist2.Replace(&falpha2);
936	//if (!fIsMonteCarlo)
937	// tlist2.Replace(&fill9b);
938	if (!fIsWobble \|\| !fNameHist.IsNull())
939	tlist2.Replace(&ffs2);
940	if (fIsWobble)
941	{
942	tlist2.AddToListAfter(&fillphi, &falpha2);
943	if (!fNameHist.IsNull())
944	tlist2.RemoveFromList(&ffs);
945	}
946
947	if (!set.IsMonteCarlo())
948	tlist.AddToList(&fillvs, "EffectiveOnTime");
949	else
950	tlist2.AddToListBefore(&fillvs, &scalc);
951	}
952
953	// by setting it here it is distributed to all consecutive tasks
954	tlist.SetAccelerator(MTask::kAccDontReset\|MTask::kAccDontTime);
955
956	par.SetVal(1);
957
958	// Execute first analysis
959	if (!evtloop.Eventloop(fMaxEvents))
960	{
961	*fLog << err << GetDescriptor() << ": Processing of on-sequences failed." << endl;
962	return -4;
963	}
964
965	if (write0)
966	delete write0;
967	if (write1)
968	delete write1;
969
970	// FIXME: Perform fit and plot energy dependant alpha plots
971	// and fit result to new tabs!
972	if (!WriteResult(plist, set))
973	return kFALSE;
974
975	*fLog << all << GetDescriptor() << ": Done." << endl;
976	*fLog << endl << endl;
977
978	return kTRUE;
979	}

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

Download in other formats: