Context Navigation

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

Visit:

Last change on this file since 8405 was 8383, checked in by tbretz, 18 years ago
* empty log message *
File size: 30.3 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)
648	return kFALSE;
649
650	const TString path(Form("%s/", fPathOut.Data()));
651	TString fname0(path);
652	TString fname1(path);
653	fname0 += fNameSummary.IsNull() ? (TString) Form("ganymed%08d-summary.root", set.GetNumAnalysis()) : fNameSummary;
654	fname1 += fNameResult.IsNull() ? (TString) Form("ganymed%08d.root", set.GetNumAnalysis()) : fNameResult;
655
656	MWriteRootFile *write0 = CanStoreSummary() ? new MWriteRootFile(fPathOut.IsNull()?0:fname0.Data(), fOverwrite?"RECREATE":"NEW") : 0;
657	MWriteRootFile *write1 = CanStoreResult() ? new MWriteRootFile(fPathOut.IsNull()?0:fname1.Data(), fOverwrite?"RECREATE":"NEW") : 0;
658	SetupWriter(write0, "WriteAfterCut0");
659	SetupWriter(write1, "WriteAfterCut3");
660
661	MTaskEnv taskenv2("CalcHadronness");
662	taskenv2.SetDefault(fCalcHadronness);
663
664	MTaskEnv taskenv3("CalcDisp");
665	taskenv3.SetDefault(fCalcDisp);
666
667	MParameterCalc setevtnum("MRawEvtHeader.fDAQEvtNumber", "SetEvtNumber");
668	setevtnum.SetNameParameter("EvtNumber");
669
670	MParameterCalc setrunnum("MRawRunHeader.fRunNumber", "SetRunNumber");
671	setrunnum.SetNameParameter("RunNumber");
672
673	// MFillH fill0a("OnPos [MHSrcPosCam]", "MSrcPosCam", "FillSrcPosCam");
674	MFillH fill1a("MHHillasOffPre [MHHillas]", "MHillas", "FillHillasPre");
675	MFillH fill2a("MHHillasOffPost [MHHillas]", "MHillas", "FillHillasPost");
676	MFillH fill3a("MHVsSizeOffPost [MHVsSize]", "MHillasSrc", "FillVsSizePost");
677	MFillH fill4a("MHHilExtOffPost [MHHillasExt]", "MHillasSrc", "FillHilExtPost");
678	MFillH fill5a("MHHilSrcOffPost [MHHillasSrc]", "MHillasSrc", "FillHilSrcPost");
679	MFillH fill6a("MHImgParOffPost [MHImagePar]", "MImagePar", "FillImgParPost");
680	MFillH fill7a("MHNewParOffPost [MHNewImagePar]", "MNewImagePar", "FillNewParPost");
681	//MFillH fill9a("MHEffOffTime [MHEffectiveOnTime]", "MTime", "FillEffOnTime");
682	fill1a.SetNameTab("PreCut");
683	fill2a.SetNameTab("PostCut");
684	fill3a.SetNameTab("VsSize");
685	fill4a.SetNameTab("HilExt");
686	fill5a.SetNameTab("HilSrc");
687	fill6a.SetNameTab("ImgPar");
688	fill7a.SetNameTab("NewPar");
689	//fill9a.SetNameTab("EffOffT");
690
691	//MFDataMember fbin("Bin.fVal", '>', 0);
692	//fill9a.SetFilter(&fbin);
693
694	MPrint print2("MEffectiveOnTime");
695	print2.EnableSkip();
696
697	// How to get source position from off- and on-data?
698	MSrcPosCalc scalc;
699	scalc.SetMode(fIsWobble?MSrcPosCalc::kWobble:MSrcPosCalc::kOffData); /********************/
700
701	MSrcPosCorrect scor;
702
703	MHillasCalc hcalc;
704	MHillasCalc hcalc2("MHillasCalcAnti");
705	hcalc.SetFlags(MHillasCalc::kCalcHillasSrc);
706	hcalc2.SetFlags(MHillasCalc::kCalcHillasSrc);
707	hcalc2.SetNameHillasSrc("MHillasSrcAnti");
708	hcalc2.SetNameSrcPosCam("MSrcPosAnti");
709
710	MSrcPosRndm srcrndm;
711
712	MTaskList tlist2;
713	tlist2.AddToList(&scalc);
714	tlist2.AddToList(&scor);
715	tlist2.AddToList(&srcrndm);
716	//if (fRndmSrcPos && !fIsWobble)
717	// tlist2.AddToList(&fill0a);
718	tlist2.AddToList(&hcalc);
719	if (fIsWobble)
720	tlist2.AddToList(&hcalc2);
721	//tlist2.AddToList(&taskenv1);
722	tlist2.AddToList(&cont0);
723	tlist2.AddToList(&taskenv2);
724	tlist2.AddToList(&taskenv3);
725	tlist2.AddToList(&setrunnum);
726	tlist2.AddToList(&setevtnum);
727	if (write0)
728	tlist2.AddToList(write0);
729	if (!fWriteOnly)
730	tlist2.AddToList(&fill1a);
731	tlist2.AddToList(&cont1);
732	if (!fWriteOnly && (!fIsWobble \|\| !fNameHistFS.IsNull()))
733	tlist2.AddToList(&ffs);
734	tlist2.AddToList(&cont2);
735	if (!fWriteOnly)
736	{
737	tlist2.AddToList(&fill2a);
738	if (fFullDisplay)
739	{
740	tlist2.AddToList(&fill3a);
741	tlist2.AddToList(&fill4a);
742	tlist2.AddToList(&fill5a);
743	tlist2.AddToList(&fill6a);
744	tlist2.AddToList(&fill7a);
745	}
746	}
747	if (!fWriteOnly)
748	{
749	tlist2.AddToList(&falpha);
750	/* if (!fIsMonteCarlo)
751	{
752	tlist2.AddToList(&fbin);
753	tlist2.AddToList(&fill9a);
754	} */
755	}
756	tlist2.AddToList(&cont3);
757	if (write1)
758	tlist2.AddToList(write1);
759
760	MPointingDevCalc devcalc;
761
762	tlist.AddToList(&readoff);
763	if (gLog.GetDebugLevel()>4)
764	tlist.AddToList(&print2, "EffectiveOnTime");
765	tlist.AddToList(&devcalc, "Starguider");
766	tlist.AddToList(&tlist2, "Events");
767
768	// by setting it here it is distributed to all consecutive tasks
769	tlist.SetAccelerator(MTask::kAccDontReset\|MTask::kAccDontTime);
770
771	par.SetVal(0);
772
773	// Create and setup the eventloop
774	MEvtLoop evtloop(fName);
775	evtloop.SetParList(&plist);
776	evtloop.SetDisplay(fDisplay);
777	evtloop.SetLogStream(fLog);
778	if (!SetupEnv(evtloop))
779	return kFALSE;
780
781	TObjArray cont;
782	cont.Add(&cont0);
783	cont.Add(&cont1);
784	cont.Add(&cont2);
785	cont.Add(&cont3);
786	//if (taskenv1.GetTask())
787	// cont.Add(taskenv1.GetTask());
788	if (taskenv2.GetTask())
789	cont.Add(taskenv2.GetTask());
790	if (taskenv3.GetTask())
791	cont.Add(taskenv3.GetTask());
792
793	if (!WriteTasks(set.GetNumAnalysis(), cont))
794	return kFALSE;
795
796	if (set.HasOffSequences() \|\| fIsWobble)
797	{
798	// Execute first analysis
799	if (!evtloop.Eventloop(fMaxEvents))
800	{
801	*fLog << err << GetDescriptor() << ": Processing of off-sequences failed." << endl;
802	return -3;
803	}
804
805	if (!evtloop.GetDisplay())
806	{
807	*fLog << err << GetDescriptor() << ": Execution stopped by user." << endl;
808	return kFALSE;
809	}
810
811	//plist.FindObject("MTimeEffectiveOnTime")->Clear();
812	}
813	else
814	{
815	// This is the simplest way to remove the two object from the parlist
816	delete halphaoff;
817	delete hfsoff;
818	}
819
820	// ------------- Loop On Data --------------------
821	MReadReports readondata;
822	readondata.AddTree("Events", "MTime.", MReadReports::kMaster);
823	readondata.AddTree("Drive", MReadReports::kRequired);
824	readondata.AddTree("Starguider", MReadReports::kRequired);
825	readondata.AddTree("EffectiveOnTime");
826
827	MReadMarsFile readonmc("Events");
828	readonmc.DisableAutoScheme();
829
830	MRead &readon = fIsMonteCarlo ? (MRead&)readonmc : (MRead&)readondata;
831	if (!set.AddFilesOn(readon))
832	return kFALSE;
833
834	scalc.SetMode(MSrcPosCalc::kDefault);
835
836	MFillH fill1b("MHHillasOnPre [MHHillas]", "MHillas", "FillHillasPre");
837	MFillH fill2b("MHHillasOnPost [MHHillas]", "MHillas", "FillHillasPost");
838	MFillH fill3b("MHVsSizeOnPost [MHVsSize]", "MHillasSrc", "FillVsSizePost");
839	MFillH fill4b("MHHilExtOnPost [MHHillasExt]", "MHillasSrc", "FillHilExtPost");
840	MFillH fill5b("MHHilSrcOnPost [MHHillasSrc]", "MHillasSrc", "FillHilSrcPost");
841	MFillH fill6b("MHImgParOnPost [MHImagePar]", "MImagePar", "FillImgParPost");
842	MFillH fill7b("MHNewParOnPost [MHNewImagePar]", "MNewImagePar", "FillNewParPost");
843	//MFillH fill9b("MHEffOnTime [MHEffectiveOnTime]", "MTime", "FillEffOnTime");
844	fill1b.SetNameTab("PreCut");
845	fill2b.SetNameTab("PostCut");
846	fill3b.SetNameTab("VsSize");
847	fill4b.SetNameTab("HilExt");
848	fill5b.SetNameTab("HilSrc");
849	fill6b.SetNameTab("ImgPar");
850	fill7b.SetNameTab("NewPar");
851	//fill9b.SetNameTab("EffOnT");
852	fill1b.SetDrawOption(set.HasOffSequences()\|\|fIsWobble?"same":"");
853	fill2b.SetDrawOption(set.HasOffSequences()\|\|fIsWobble?"same":"");
854	fill3b.SetDrawOption(set.HasOffSequences()\|\|fIsWobble?"same":"");
855	fill4b.SetDrawOption(set.HasOffSequences()\|\|fIsWobble?"same":"");
856	fill5b.SetDrawOption(set.HasOffSequences()\|\|fIsWobble?"same":"");
857	fill6b.SetDrawOption(set.HasOffSequences()\|\|fIsWobble?"same":"");
858	fill7b.SetDrawOption(set.HasOffSequences()\|\|fIsWobble?"same":"");
859
860	//fill9b.SetFilter(&fbin);
861
862	/*
863	MHVsTime hvs("MEffectiveOnTime.fVal");
864	hvs.SetTitle("Effective On-Time vs. Time;;T_{on}");
865	MFillH fillvs(&hvs, "MTimeEffectiveOnTime", "FillOnTime");
866	fillvs.SetNameTab("OnTime");
867	*/
868	MH3 hvs("MPointingPos.fZd");
869	hvs.SetName("Theta");
870	hvs.SetTitle("Effective On-Time vs. Zenith Angle;\\Theta [\\circ];T_{on} [s]");
871	MFillH fillvs(&hvs, "", "FillOnTime");
872	if (!fIsMonteCarlo)
873	fillvs.SetWeight("MEffectiveOnTime");
874	fillvs.SetNameTab("OnTime");
875
876	/*
877	MParameterD weight;
878	weight.SetVal(-1);
879	fill2a.SetWeight(&weight);
880	fill3a.SetWeight(&weight);
881	fill4a.SetWeight(&weight);
882	fill5a.SetWeight(&weight);
883	fill6a.SetWeight(&weight);
884	fill7a.SetWeight(&weight);
885	if (fSubstraction)
886	{
887	fill2a.SetNameTab("PostCut-");
888	fill3a.SetNameTab("VsSize-");
889	fill4a.SetNameTab("HilExt-");
890	fill5a.SetNameTab("HilSrc-");
891	fill6a.SetNameTab("ImgPar-");
892	fill7a.SetNameTab("NewPar-");
893	}
894	*/
895	MHAlpha *halphaon=CreateNewHist(plist);
896	MFillH falpha2(halphaon, "", "FillHist");
897	MH *hfs=CreateNewHistFS(plist);
898	MFillH ffs2(hfs, "MHillas", "FillFS");
899	MFillH fillphi("MHPhi", "", "FillPhi");
900	fillphi.SetDrawOption("anticut");
901
902	tlist.Replace(&readon);
903	if (fRndmSrcPos && !fIsWobble)
904	{
905	// tlist2.RemoveFromList(&fill0a);
906	tlist2.RemoveFromList(&srcrndm);
907	}
908
909	MFillH fillsrc(&hsrcpos, "MSrcPosCam", "FillSrcPosCam");
910	fillsrc.SetNameTab("SrcPos");
911	if (fIsWobble)
912	tlist2.AddToListBefore(&fillsrc, &hcalc);
913
914	if (!fWriteOnly)
915	{
916	tlist2.Replace(&fill1b);
917	/* if (fIsWobble)
918	{
919	tlist2.AddToListAfter(&fill2b, &fill2a);
920	tlist2.AddToListAfter(&fill3b, &fill3a);
921	}
922	else
923	*/
924	tlist2.Replace(&fill2b);
925	if (fFullDisplay)
926	{
927	tlist2.Replace(&fill3b);
928	tlist2.Replace(&fill4b);
929	tlist2.Replace(&fill5b);
930	tlist2.Replace(&fill6b);
931	tlist2.Replace(&fill7b);
932	}
933	tlist2.Replace(&falpha2);
934	//if (!fIsMonteCarlo)
935	// tlist2.Replace(&fill9b);
936	if (!fIsWobble \|\| !fNameHist.IsNull())
937	tlist2.Replace(&ffs2);
938	if (fIsWobble)
939	{
940	tlist2.AddToListAfter(&fillphi, &falpha2);
941	if (!fNameHist.IsNull())
942	tlist2.RemoveFromList(&ffs);
943	}
944
945	if (!fIsMonteCarlo)
946	tlist.AddToList(&fillvs, "EffectiveOnTime");
947	else
948	tlist2.AddToListBefore(&fillvs, &scalc);
949	}
950
951	// by setting it here it is distributed to all consecutive tasks
952	tlist.SetAccelerator(MTask::kAccDontReset\|MTask::kAccDontTime);
953
954	par.SetVal(1);
955
956	// Execute first analysis
957	if (!evtloop.Eventloop(fMaxEvents))
958	{
959	*fLog << err << GetDescriptor() << ": Processing of on-sequences failed." << endl;
960	return -4;
961	}
962
963	if (write0)
964	delete write0;
965	if (write1)
966	delete write1;
967
968	// FIXME: Perform fit and plot energy dependant alpha plots
969	// and fit result to new tabs!
970	if (!WriteResult(plist, set))
971	return kFALSE;
972
973	*fLog << all << GetDescriptor() << ": Done." << endl;
974	*fLog << endl << endl;
975
976	return kTRUE;
977	}

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

Download in other formats: