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

Last change on this file since 6901 was 6892, checked in by tbretz, 20 years ago
*** empty log message ***
File size: 15.9 KB
Line 
1/* ======================================================================== *\
2!
3! *
4! * This file is part of MARS, the MAGIC Analysis and Reconstruction
5! * Software. It is distributed to you in the hope that it can be a useful
6! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
7! * It is distributed WITHOUT ANY WARRANTY.
8! *
9! * Permission to use, copy, modify and distribute this software and its
10! * documentation for any purpose is hereby granted without fee,
11! * provided that the above copyright notice appear in all copies and
12! * that both that copyright notice and this permission notice appear
13! * in supporting documentation. It is provided "as is" without express
14! * or implied warranty.
15! *
16!
17!
18! Author(s): Thomas Bretz, 1/2005 <mailto:tbretz@astro.uni-wuerzburg.de>
19!
20! Copyright: MAGIC Software Development, 2000-2005
21!
22!
23\* ======================================================================== */
24
25/////////////////////////////////////////////////////////////////////////////
26//
27// MJCut
28//
29// FIXME: Preparation for wobble mode missing
30//
31/////////////////////////////////////////////////////////////////////////////
32#include "MJCut.h"
33
34#include <TEnv.h>
35#include <TFile.h>
36
37#include "MLog.h"
38#include "MLogManip.h"
39
40#include "MParList.h"
41#include "MTaskList.h"
42#include "MEvtLoop.h"
43
44#include "MStatusDisplay.h"
45
46#include "MReadReports.h"
47#include "MPrint.h"
48#include "MContinue.h"
49#include "MEnergyEstimate.h"
50#include "MTaskEnv.h"
51#include "MSrcPosCalc.h"
52#include "MHillasCalc.h"
53#include "MFillH.h"
54#include "MWriteRootFile.h"
55
56#include "../mhflux/MAlphaFitter.h"
57#include "MBinning.h"
58#include "MDataSet.h"
59#include "MParameters.h"
60#include "MPointingPos.h"
61#include "MObservatory.h"
62
63ClassImp(MJCut);
64
65using namespace std;
66
67// --------------------------------------------------------------------------
68//
69// Default constructor.
70//
71MJCut::MJCut(const char *name, const char *title)
72 : fStoreSummary(kFALSE), fStoreResult(kFALSE), fWriteOnly(kFALSE),
73 fIsWobble(kFALSE), fEstimateEnergy(0), fCalcHadronness(0)
74{
75 fName = name ? name : "MJCut";
76 fTitle = title ? title : "Standard program to perform g/h-seperation cuts";
77}
78
79MJCut::~MJCut()
80{
81 if (fEstimateEnergy)
82 delete fEstimateEnergy;
83 if (fCalcHadronness)
84 delete fCalcHadronness;
85}
86
87// --------------------------------------------------------------------------
88//
89// Set the name of the summary file (events after cut0)
90// If you give a name the storage of this file is enabled implicitly.
91// If you give no filename the storage is neither enabled nor disabled,
92// but the storage file name is reset.
93// If no filename is set the default filename is used.
94// You can explicitly enable or disable the storage using EnableStoreOf*()
95// The default argument is no filename.
96//
97void MJCut::SetNameSummaryFile(const char *name)
98{
99 fNameSummary=name;
100 if (!fNameSummary.IsNull())
101 EnableStorageOfSummary();
102}
103
104// --------------------------------------------------------------------------
105//
106// Set the name of the summary file (events after cut3)
107// If you give a name the storage of this file is enabled implicitly.
108// If you give no filename the storage is neither enabled nor disabled,
109// but the storage file name is reset.
110// If no filename is set the default filename is used.
111// You can explicitly enable or disable the storage using EnableStoreOf*()
112// The default argument is no filename.
113//
114void MJCut::SetNameResultFile(const char *name)
115{
116 fNameResult=name;
117 if (!fNameResult.IsNull())
118 EnableStorageOfResult();
119}
120
121// --------------------------------------------------------------------------
122//
123// Setup a task estimating the energy. The given task is cloned.
124//
125void MJCut::SetEnergyEstimator(const MTask *task)
126{
127 if (fEstimateEnergy)
128 delete fEstimateEnergy;
129 fEstimateEnergy = task ? (MTask*)task->Clone() : 0;
130}
131
132// --------------------------------------------------------------------------
133//
134// Setup a task calculating the hadronness. The given task is cloned.
135//
136void MJCut::SetHadronnessCalculator(const MTask *task)
137{
138 if (fCalcHadronness)
139 delete fCalcHadronness;
140 fCalcHadronness = task ? (MTask*)task->Clone() : 0;
141}
142
143// --------------------------------------------------------------------------
144//
145// return fOutputPath+"/ganymed%08d.root", num
146//
147TString MJCut::GetOutputFile(UInt_t num) const
148{
149 TString p(fPathOut);
150 p += "/";
151 p += fNameOutput.IsNull() ? Form("ganymed%08d.root", num) : fNameOutput.Data();
152 return p;
153}
154
155/*
156Bool_t MJCut::ReadTasks(const char *fname, MTask* &env1, MTask* &env2) const
157{
158 // const TString fname = Form("%s/calib%08d.root", fPathIn.Data(), fSequence.GetSequence());
159
160 *fLog << inf << "Reading from file: " << fname << endl;
161
162 TFile file(fname, "READ");
163 if (!file.IsOpen())
164 {
165 *fLog << err << dbginf << "ERROR - Could not open file " << fname << endl;
166 return kFALSE;
167 }
168
169 TObject *o = file.Get("EstimateEnergy");
170 if (o && !o->InheritsFrom(MTask::Class()))
171 {
172 *fLog << err << dbginf << "ERROR - EstimateEnergy read from " << fname << " doesn't inherit from MTask!" << endl;
173 return kFALSE;
174 }
175 env1 = o ? (MTask*)o->Clone() : NULL;
176
177 o = file.Get("CalcHadronness");
178 if (o && !o->InheritsFrom(MTask::Class()))
179 {
180 *fLog << err << dbginf << "ERROR - CalcHadronness read from " << fname << " doesn't inherit from MTask!" << endl;
181 return kFALSE;
182 }
183 env2 = o ? (MTask*)o->Clone() : NULL;
184
185 return kTRUE;
186}
187*/
188
189// --------------------------------------------------------------------------
190//
191// Write the tasks in cont to the file corresponding to analysis number num,
192// see GetOutputFile()
193//
194Bool_t MJCut::WriteTasks(UInt_t num, TObjArray &cont) const
195{
196 if (fPathOut.IsNull())
197 {
198 *fLog << inf << "No output path specified via SetPathOut - no output written." << endl;
199 return kTRUE;
200 }
201
202 const TString oname(GetOutputFile(num));
203
204 *fLog << inf << "Writing to file: " << oname << endl;
205
206 TFile file(oname, fOverwrite?"RECREATE":"NEW", "File created by MJCut", 9);
207 if (!file.IsOpen())
208 {
209 *fLog << err << "ERROR - Couldn't open file " << oname << " for writing..." << endl;
210 return kFALSE;
211 }
212
213 return WriteContainer(cont);
214}
215
216// --------------------------------------------------------------------------
217//
218// Write the result plots and other results to the file corresponding to
219// analysis number num, see GetOutputFile()
220//
221Bool_t MJCut::WriteResult(UInt_t num) const
222{
223 TObjArray arr;
224 return WriteContainer(arr, GetOutputFile(num), "UPDATE");
225}
226
227// --------------------------------------------------------------------------
228//
229// MJCut allows to setup several option by a resource file:
230// MJCut.WriteSummary: yes, no
231// MJCut.SummaryFile: filename
232// MJCut.WriteResult: yes, no
233// MJCut.ResultFile: filename
234//
235Bool_t MJCut::CheckEnvLocal()
236{
237 const TString f0(GetEnv("SummaryFile", ""));
238 const TString f1(GetEnv("ResultFile", ""));
239 if (!f0.IsNull())
240 SetNameSummaryFile(f0);
241 if (!f1.IsNull())
242 SetNameResultFile(f1);
243
244 EnableStorageOfSummary(GetEnv("SummaryFile", fStoreSummary));
245 EnableStorageOfResult(GetEnv("ResultFile", fStoreResult));
246 EnableWobbleMode(GetEnv("WobbleMode", fIsWobble));
247
248 return kTRUE;
249}
250
251// --------------------------------------------------------------------------
252//
253// Setup write to write:
254// container tree optional?
255// -------------- ---------- -----------
256// "MHillas" to "Events"
257// "MHillasSrc" to "Events"
258// "MHadronness" to "Events" yes
259// "MEnergyEst" to "Events" yes
260// "DataType" to "Events"
261//
262void MJCut::SetupWriter(MWriteRootFile &write, const char *name) const
263{
264 write.SetName(name);
265 write.AddContainer("MHillas", "Events");
266 write.AddContainer("MHillasSrc", "Events");
267 write.AddContainer("MHillasSrcAnti", "Events", kFALSE);
268 write.AddContainer("MHadronness", "Events", kFALSE);
269 write.AddContainer("MEnergyEst", "Events", kFALSE);
270 write.AddContainer("DataType", "Events");
271
272 // Should not be the default: Either as option, or as
273 // setup from resource file
274 // write.AddContainer("MHillasExt", "Events");
275 // write.AddContainer("MImagePar", "Events");
276 // write.AddContainer("MNewImagePar", "Events");
277}
278
279Bool_t MJCut::ProcessFile(const MDataSet &set)
280{
281 if (!set.IsValid())
282 {
283 *fLog << err << "ERROR - DataSet invalid!" << endl;
284 return kFALSE;
285 }
286
287 CheckEnv();
288
289 // --------------------------------------------------------------------------------
290
291 *fLog << inf;
292 fLog->Separator(GetDescriptor());
293 *fLog << "Perform cuts for data set " << set.GetName() << endl;
294 *fLog << endl;
295
296 // --------------------------------------------------------------------------------
297
298 // Setup Parlist
299 MParList plist;
300 plist.AddToList(this); // take care of fDisplay!
301
302 MParameterI par("DataType");
303 plist.AddToList(&par);
304
305 // Setup Tasklist
306 MTaskList tlist;
307 plist.AddToList(&tlist);
308
309 // La Palma Magic1
310 MObservatory obs;
311 plist.AddToList(&obs);
312
313 // Possible source position (eg. Wobble Mode)
314 MPointingPos source("MSourcePos");
315 if (set.GetSourcePos(source))
316 {
317 plist.AddToList(&source);
318 *fLog << inf << "Using Source Position: " << source.GetTitle() << endl;
319 }
320 else
321 *fLog << inf << "No source position applied..." << endl;
322
323 // Initialize default binnings
324 MBinning bins1(18, 0, 90, "BinningAlpha", "lin");
325 MBinning bins2(15, 10, 1e6 , "BinningEnergyEst", "log");
326 MBinning bins3(50, 0, 60, "BinningTheta", "cos");
327 MBinning bins4("BinningFalseSource");
328 plist.AddToList(&bins1);
329 plist.AddToList(&bins2);
330 plist.AddToList(&bins3);
331 plist.AddToList(&bins4);
332
333 // --------------------------------------------------------------------------------
334
335 // Setup fitter and histograms
336 MAlphaFitter fit;
337 plist.AddToList(&fit);
338 if (fIsWobble)
339 fit.SetScaleMode(MAlphaFitter::kNone);
340
341 MFillH falpha("MHAlphaOff [MHAlpha]", "MHillasSrc", "FillAlpha");
342 MFillH ffs("MHFalseSourceOff [MHFalseSource]", "MHillas", "FillFS");
343
344 // FIXME: If fPathIn read cuts and energy estimator from file!
345 MContinue cont0("", "Cut0");
346 MContinue cont1("", "Cut1");
347 MContinue cont2("", "Cut2");
348 MContinue cont3("", "Cut3");
349 cont0.SetAllowEmpty();
350 cont1.SetAllowEmpty();
351 cont2.SetAllowEmpty();
352 cont3.SetAllowEmpty();
353
354 // ------------- Loop Off Data --------------------
355 MReadReports readoff;
356 readoff.AddTree("Events", "MTime.", kTRUE);
357 readoff.AddTree("Drive");
358 readoff.AddTree("EffectiveOnTime");
359 if (fIsWobble)
360 set.AddFilesOn(readoff);
361 else
362 set.AddFilesOff(readoff);
363
364 const TString path(Form("%s/", fPathOut.Data()));
365 TString fname0(path);
366 TString fname1(path);
367 fname0 += fNameSummary.IsNull() ? (TString) Form("ganymed%08d-summary.root", set.GetNumAnalysis()) : fNameSummary;
368 fname1 += fNameResult.IsNull() ? (TString) Form("ganymed%08d-result.root", set.GetNumAnalysis()) : fNameResult;
369
370 MWriteRootFile write0(fPathOut.IsNull()?0:fname0.Data(), fOverwrite?"RECREATE":"NEW");
371 MWriteRootFile write1(fPathOut.IsNull()?0:fname1.Data(), fOverwrite?"RECREATE":"NEW");
372 if (CanStoreSummary())
373 SetupWriter(write0, "WriteAfterCut0");
374 if (CanStoreSummary())
375 SetupWriter(write1, "WriteAfterCut3");
376
377
378 MEnergyEstimate est;
379
380 MTaskEnv taskenv1("EstimateEnergy");
381 taskenv1.SetDefault(fEstimateEnergy ? fEstimateEnergy : &est);
382
383 MTaskEnv taskenv2("CalcHadronness");
384 taskenv2.SetDefault(fCalcHadronness);
385
386 MFillH fill1a("MHHillasOffPre [MHHillas]", "MHillas", "FillHillasPre");
387 MFillH fill2a("MHHillasOffPost [MHHillas]", "MHillas", "FillHillasPost");
388 MFillH fill3a("MHVsSizeOffPost [MHVsSize]", "MHillasSrc", "FillVsSizePost");
389 fill1a.SetNameTab("PreCut");
390 fill2a.SetNameTab("PostCut");
391 fill3a.SetNameTab("VsSize");
392
393 MPrint print2("MEffectiveOnTime");
394
395 // How to get source position from off- and on-data?
396 MSrcPosCalc scalc;
397 if (fIsWobble)
398 scalc.SetWobbleMode(); /********************/
399 MHillasCalc hcalc;
400 MHillasCalc hcalc2("MHillasCalcAnti");
401 hcalc.SetFlags(MHillasCalc::kCalcHillasSrc);
402 hcalc2.SetFlags(MHillasCalc::kCalcHillasSrc);
403 hcalc2.SetNameHillasSrc("MHillasSrcAnti");
404 hcalc2.SetNameSrcPosCam("MSrcPosAnti");
405
406 MTaskList tlist2;
407 tlist2.AddToList(&scalc);
408 tlist2.AddToList(&hcalc);
409 if (fIsWobble)
410 tlist2.AddToList(&hcalc2);
411 tlist2.AddToList(&taskenv1);
412 tlist2.AddToList(&taskenv2);
413 tlist2.AddToList(&cont0);
414 if (CanStoreSummary())
415 tlist2.AddToList(&write0);
416 if (!fWriteOnly)
417 tlist2.AddToList(&fill1a);
418 tlist2.AddToList(&cont1);
419 if (!fWriteOnly)
420 tlist2.AddToList(&ffs);
421 tlist2.AddToList(&cont2);
422 if (!fWriteOnly)
423 {
424 tlist2.AddToList(&fill2a);
425 tlist2.AddToList(&fill3a);
426 }
427 if (!fWriteOnly)
428 tlist2.AddToList(&falpha);
429 tlist2.AddToList(&cont3);
430 if (CanStoreResult())
431 tlist2.AddToList(&write1);
432
433 tlist.AddToList(&readoff);
434 tlist.AddToList(&print2, "EffectiveOnTime");
435 tlist.AddToList(&tlist2, "Events");
436
437 par.SetVal(0);
438
439 // Create and setup the eventloop
440 MEvtLoop evtloop(fName);
441 evtloop.SetParList(&plist);
442 evtloop.SetDisplay(fDisplay);
443 evtloop.SetLogStream(fLog);
444 if (!SetupEnv(evtloop))
445 return kFALSE;
446
447 if (set.HasOffSequences() || fIsWobble)
448 {
449 // Execute first analysis
450 if (!evtloop.Eventloop(fMaxEvents))
451 {
452 *fLog << err << GetDescriptor() << ": Processing of off-sequences failed." << endl;
453 return kFALSE;
454 }
455
456 tlist.PrintStatistics();
457
458 if (!evtloop.GetDisplay())
459 {
460 *fLog << err << GetDescriptor() << ": Execution stopped by user." << endl;
461 return kFALSE;
462 }
463 }
464
465 // ------------- Loop On Data --------------------
466 MReadReports readon;
467 readon.AddTree("Events", "MTime.", kTRUE);
468 readon.AddTree("Drive");
469 readon.AddTree("EffectiveOnTime");
470 set.AddFilesOn(readon);
471
472 if (fIsWobble)
473 scalc.SetWobbleMode(kFALSE); /********************/
474
475 MFillH fill1b("MHHillasOnPre [MHHillas]", "MHillas", "FillHillasPre");
476 MFillH fill2b("MHHillasOnPost [MHHillas]", "MHillas", "FillHillasPost");
477 MFillH fill3b("MHVsSizeOnPost [MHVsSize]", "MHillasSrc", "FillVsSizePost");
478 fill1b.SetNameTab("PreCut");
479 fill2b.SetNameTab("PostCut");
480 fill3b.SetNameTab("VsSize");
481 fill1b.SetDrawOption(set.HasOffSequences()||fIsWobble?"same":"");
482 fill2b.SetDrawOption(set.HasOffSequences()||fIsWobble?"same":"");
483 fill3b.SetDrawOption(set.HasOffSequences()||fIsWobble?"same":"");
484
485 MFillH falpha2("MHAlpha", "MHillasSrc", "FillAlpha");
486 MFillH ffs2("MHFalseSource", "MHillas", "FillFS");
487
488 tlist.Replace(&readon);
489 if (!fWriteOnly)
490 {
491 tlist2.Replace(&fill1b);
492 tlist2.Replace(&fill2b);
493 tlist2.Replace(&fill3b);
494 tlist2.Replace(&falpha2);
495 tlist2.Replace(&ffs2);
496 }
497
498 par.SetVal(1);
499
500 TObjArray cont;
501 cont.Add(&cont0);
502 cont.Add(&cont1);
503 cont.Add(&cont2);
504 cont.Add(&cont3);
505 if (taskenv1.GetTask())
506 cont.Add(taskenv1.GetTask());
507 if (taskenv2.GetTask())
508 cont.Add(taskenv2.GetTask());
509
510 if (!WriteTasks(set.GetNumAnalysis(), cont))
511 return kFALSE;
512
513 // Execute first analysis
514 if (!evtloop.Eventloop(fMaxEvents))
515 {
516 *fLog << err << GetDescriptor() << ": Processing of on-sequences failed." << endl;
517 return kFALSE;
518 }
519
520 tlist.PrintStatistics();
521
522 // FIXME: Perform fit and plot energy dependant alpha plots
523 // and fit result to new tabs!
524 if (!WriteResult(set.GetNumAnalysis()))
525 return kFALSE;
526
527 *fLog << all << GetDescriptor() << ": Done." << endl;
528 *fLog << endl << endl;
529
530 return kTRUE;
531}
Note: See TracBrowser for help on using the repository browser.