Context Navigation

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

Visit:

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

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

Download in other formats: