Context Navigation

MJCut.cc@ 19689

Visit:

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

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

Context Navigation

source: trunk/Mars/mjobs/MJCut.cc@ 19689

Download in other formats: