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

Last change on this file since 6877 was 6874, checked in by tbretz, 20 years ago
*** empty log message ***
File size: 15.4 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("MHadronness", "Events", kFALSE);
268 write.AddContainer("MEnergyEst", "Events", kFALSE);
269 write.AddContainer("DataType", "Events");
270
271 // Should not be the default: Either as option, or as
272 // setup from resource file
273 // write.AddContainer("MHillasExt", "Events");
274 // write.AddContainer("MImagePar", "Events");
275 // write.AddContainer("MNewImagePar", "Events");
276}
277
278Bool_t MJCut::ProcessFile(const MDataSet &set)
279{
280 if (!set.IsValid())
281 {
282 *fLog << err << "ERROR - DataSet invalid!" << endl;
283 return kFALSE;
284 }
285
286 CheckEnv();
287
288 // --------------------------------------------------------------------------------
289
290 *fLog << inf;
291 fLog->Separator(GetDescriptor());
292 *fLog << "Perform cuts for data set " << set.GetName() << endl;
293 *fLog << endl;
294
295 // --------------------------------------------------------------------------------
296
297 // Setup Parlist
298 MParList plist;
299 plist.AddToList(this); // take care of fDisplay!
300
301 MParameterI par("DataType");
302 plist.AddToList(&par);
303
304 // Setup Tasklist
305 MTaskList tlist;
306 plist.AddToList(&tlist);
307
308 // La Palma Magic1
309 MObservatory obs;
310 plist.AddToList(&obs);
311
312 // Possible source position (eg. Wobble Mode)
313 MPointingPos source("MSourcePos");
314 if (set.GetSourcePos(source))
315 {
316 plist.AddToList(&source);
317 *fLog << inf << "Using Source Position: " << source.GetTitle() << endl;
318 }
319 else
320 *fLog << inf << "No source position applied..." << endl;
321
322 // Initialize default binnings
323 MBinning bins1(18, 0, 90, "BinningAlpha", "lin");
324 MBinning bins2(25, 10, 1e6 , "BinningEnergyEst", "log");
325 MBinning bins3(50, 0, 60, "BinningTheta", "cos");
326 MBinning bins4("BinningFalseSource");
327 plist.AddToList(&bins1);
328 plist.AddToList(&bins2);
329 plist.AddToList(&bins3);
330 plist.AddToList(&bins4);
331
332 // --------------------------------------------------------------------------------
333
334 // Setup fitter and histograms
335 MAlphaFitter fit;
336 plist.AddToList(&fit);
337 if (fIsWobble)
338 fit.SetScaleMode(MAlphaFitter::kNone);
339
340 MFillH falpha("MHAlphaOff [MHAlpha]", "MHillasSrc", "FillAlpha");
341 MFillH ffs("MHFalseSourceOff [MHFalseSource]", "MHillas", "FillFS");
342
343 // FIXME: If fPathIn read cuts and energy estimator from file!
344 MContinue cont0("", "Cut0");
345 MContinue cont1("", "Cut1");
346 MContinue cont2("", "Cut2");
347 MContinue cont3("", "Cut3");
348 cont0.SetAllowEmpty();
349 cont1.SetAllowEmpty();
350 cont2.SetAllowEmpty();
351 cont3.SetAllowEmpty();
352
353 // ------------- Loop Off Data --------------------
354 MReadReports readoff;
355 readoff.AddTree("Events", "MTime.", kTRUE);
356 readoff.AddTree("Drive");
357 readoff.AddTree("EffectiveOnTime");
358 if (fIsWobble)
359 set.AddFilesOn(readoff);
360 else
361 set.AddFilesOff(readoff);
362
363 const TString path(Form("%s/", fPathOut.Data()));
364 TString fname0(path);
365 TString fname1(path);
366 fname0 += fNameSummary.IsNull() ? (TString) Form("ganymed%08d-summary.root", set.GetNumAnalysis()) : fNameSummary;
367 fname1 += fNameResult.IsNull() ? (TString) Form("ganymed%08d-result.root", set.GetNumAnalysis()) : fNameResult;
368
369 MWriteRootFile write0(fPathOut.IsNull()?0:fname0.Data(), fOverwrite?"RECREATE":"NEW");
370 MWriteRootFile write1(fPathOut.IsNull()?0:fname1.Data(), fOverwrite?"RECREATE":"NEW");
371 if (CanStoreSummary())
372 SetupWriter(write0, "WriteAfterCut0");
373 if (CanStoreSummary())
374 SetupWriter(write1, "WriteAfterCut3");
375
376
377 MEnergyEstimate est;
378
379 MTaskEnv taskenv1("EstimateEnergy");
380 taskenv1.SetDefault(fEstimateEnergy ? fEstimateEnergy : &est);
381
382 MTaskEnv taskenv2("CalcHadronness");
383 taskenv2.SetDefault(fCalcHadronness);
384
385 MFillH fill1a("MHHillasOffPre [MHHillas]", "MHillas", "FillHillasPre");
386 MFillH fill2a("MHHillasOffPost [MHHillas]", "MHillas", "FillHillasPost");
387 fill1a.SetNameTab("PreCut");
388 fill2a.SetNameTab("PostCut");
389
390 MPrint print2("MEffectiveOnTime");
391
392 // How to get source position from off- and on-data?
393 MSrcPosCalc scalc;
394 if (fIsWobble)
395 scalc.SetWobbleMode(); /********************/
396 MHillasCalc hcalc;
397 MHillasCalc hcalc2("MHillasCalcAnti");
398 hcalc.SetFlags(MHillasCalc::kCalcHillasSrc);
399 hcalc2.SetFlags(MHillasCalc::kCalcHillasSrc);
400 hcalc2.SetNameHillasSrc("MHillasSrcAnti");
401 hcalc2.SetNameSrcPosCam("MSrcPosAnti");
402
403 MTaskList tlist2;
404 tlist2.AddToList(&scalc);
405 tlist2.AddToList(&hcalc);
406 if (fIsWobble)
407 tlist2.AddToList(&hcalc2);
408 tlist2.AddToList(&taskenv1);
409 tlist2.AddToList(&taskenv2);
410 tlist2.AddToList(&cont0);
411 if (CanStoreSummary())
412 tlist2.AddToList(&write0);
413 if (!fWriteOnly)
414 tlist2.AddToList(&fill1a);
415 tlist2.AddToList(&cont1);
416 if (!fWriteOnly)
417 tlist2.AddToList(&ffs);
418 tlist2.AddToList(&cont2);
419 if (!fWriteOnly)
420 tlist2.AddToList(&fill2a);
421 if (!fWriteOnly)
422 tlist2.AddToList(&falpha);
423 tlist2.AddToList(&cont3);
424 if (CanStoreResult())
425 tlist2.AddToList(&write1);
426
427 tlist.AddToList(&readoff);
428 tlist.AddToList(&print2, "EffectiveOnTime");
429 tlist.AddToList(&tlist2, "Events");
430
431 par.SetVal(0);
432
433 // Create and setup the eventloop
434 MEvtLoop evtloop(fName);
435 evtloop.SetParList(&plist);
436 evtloop.SetDisplay(fDisplay);
437 evtloop.SetLogStream(fLog);
438 if (!SetupEnv(evtloop))
439 return kFALSE;
440
441 if (set.HasOffSequences() || fIsWobble)
442 {
443 // Execute first analysis
444 if (!evtloop.Eventloop(fMaxEvents))
445 {
446 *fLog << err << GetDescriptor() << ": Processing of off-sequences failed." << endl;
447 return kFALSE;
448 }
449
450 tlist.PrintStatistics();
451
452 if (!evtloop.GetDisplay())
453 {
454 *fLog << err << GetDescriptor() << ": Execution stopped by user." << endl;
455 return kFALSE;
456 }
457 }
458
459 // ------------- Loop On Data --------------------
460 MReadReports readon;
461 readon.AddTree("Events", "MTime.", kTRUE);
462 readon.AddTree("Drive");
463 readon.AddTree("EffectiveOnTime");
464 set.AddFilesOn(readon);
465
466 if (fIsWobble)
467 scalc.SetWobbleMode(kFALSE); /********************/
468
469 MFillH fill1b("MHHillasOnPre [MHHillas]", "MHillas", "FillHillasPre");
470 MFillH fill2b("MHHillasOnPost [MHHillas]", "MHillas", "FillHillasPost");
471 fill1b.SetNameTab("PreCut");
472 fill2b.SetNameTab("PostCut");
473 fill1b.SetDrawOption(set.HasOffSequences()||fIsWobble?"same":"");
474 fill2b.SetDrawOption(set.HasOffSequences()||fIsWobble?"same":"");
475
476 MFillH falpha2("MHAlpha", "MHillasSrc", "FillAlpha");
477 MFillH ffs2("MHFalseSource", "MHillas", "FillFS");
478
479 tlist.Replace(&readon);
480 if (!fWriteOnly)
481 {
482 tlist2.Replace(&fill1b);
483 tlist2.Replace(&fill2b);
484 tlist2.Replace(&falpha2);
485 tlist2.Replace(&ffs2);
486 }
487
488 par.SetVal(1);
489
490 TObjArray cont;
491 cont.Add(&cont0);
492 cont.Add(&cont1);
493 cont.Add(&cont2);
494 cont.Add(&cont3);
495 if (taskenv1.GetTask())
496 cont.Add(taskenv1.GetTask());
497 if (taskenv2.GetTask())
498 cont.Add(taskenv2.GetTask());
499
500 if (!WriteTasks(set.GetNumAnalysis(), cont))
501 return kFALSE;
502
503 // Execute first analysis
504 if (!evtloop.Eventloop(fMaxEvents))
505 {
506 *fLog << err << GetDescriptor() << ": Processing of on-sequences failed." << endl;
507 return kFALSE;
508 }
509
510 tlist.PrintStatistics();
511
512 // FIXME: Perform fit and plot energy dependant alpha plots
513 // and fit result to new tabs!
514 if (!WriteResult(set.GetNumAnalysis()))
515 return kFALSE;
516
517 *fLog << all << GetDescriptor() << ": Done." << endl;
518 *fLog << endl << endl;
519
520 return kTRUE;
521}
Note: See TracBrowser for help on using the repository browser.