Context Navigation

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

Visit:

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

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

Download in other formats: