Context Navigation

source: tags/Mars-V2.1/mjobs/MJCut.cc

Visit:

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

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

Download in other formats: