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

Last change on this file since 6919 was 6917, checked in by tbretz, 21 years ago
*** empty log message ***
File size: 19.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-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 "MH3.h"
58#include "MBinning.h"
59#include "MDataSet.h"
60#include "MParameters.h"
61#include "MPointingPos.h"
62#include "MObservatory.h"
63
64ClassImp(MJCut);
65
66using namespace std;
67
68// --------------------------------------------------------------------------
69//
70// Default constructor. Set defaults for fStoreSummary, fStoreresult,
71// fWriteOnly, fIsWobble and fFullDisplay to kFALSE and initialize
72// fEstimateEnergy and fCalcHadronness with NULL.
73//
74MJCut::MJCut(const char *name, const char *title)
75 : fStoreSummary(kFALSE), fStoreResult(kFALSE), fWriteOnly(kFALSE),
76 fIsWobble(kFALSE), fFullDisplay(kFALSE), fEstimateEnergy(0),
77 fCalcHadronness(0)
78{
79 fName = name ? name : "MJCut";
80 fTitle = title ? title : "Standard program to perform g/h-separation cuts";
81}
82
83// --------------------------------------------------------------------------
84//
85// Destructor. Delete fEstimateEnergy and fCalcHadronness if != NULL
86//
87MJCut::~MJCut()
88{
89 if (fEstimateEnergy)
90 delete fEstimateEnergy;
91 if (fCalcHadronness)
92 delete fCalcHadronness;
93}
94
95// --------------------------------------------------------------------------
96//
97// Set the name of the summary file (events after cut0)
98// If you give a name the storage of this file is enabled implicitly.
99// If you give no filename the storage is neither enabled nor disabled,
100// but the storage file name is reset.
101// If no filename is set the default filename is used.
102// You can explicitly enable or disable the storage using EnableStoreOf*()
103// The default argument is no filename.
104//
105void MJCut::SetNameSummaryFile(const char *name)
106{
107 fNameSummary=name;
108 if (!fNameSummary.IsNull())
109 EnableStorageOfSummary();
110}
111
112// --------------------------------------------------------------------------
113//
114// Set the name of the summary file (events after cut3)
115// If you give a name the storage of this file is enabled implicitly.
116// If you give no filename the storage is neither enabled nor disabled,
117// but the storage file name is reset.
118// If no filename is set the default filename is used.
119// You can explicitly enable or disable the storage using EnableStoreOf*()
120// The default argument is no filename.
121//
122void MJCut::SetNameResultFile(const char *name)
123{
124 fNameResult=name;
125 if (!fNameResult.IsNull())
126 EnableStorageOfResult();
127}
128
129// --------------------------------------------------------------------------
130//
131// Setup a task estimating the energy. The given task is cloned.
132//
133void MJCut::SetEnergyEstimator(const MTask *task)
134{
135 if (fEstimateEnergy)
136 delete fEstimateEnergy;
137 fEstimateEnergy = task ? (MTask*)task->Clone() : 0;
138}
139
140// --------------------------------------------------------------------------
141//
142// Setup a task calculating the hadronness. The given task is cloned.
143//
144void MJCut::SetHadronnessCalculator(const MTask *task)
145{
146 if (fCalcHadronness)
147 delete fCalcHadronness;
148 fCalcHadronness = task ? (MTask*)task->Clone() : 0;
149}
150
151// --------------------------------------------------------------------------
152//
153// return fOutputPath+"/ganymed%08d.root", num
154//
155TString MJCut::GetOutputFile(UInt_t num) const
156{
157 TString p(fPathOut);
158 p += "/";
159 p += fNameOutput.IsNull() ? Form("ganymed%08d.root", num) : fNameOutput.Data();
160 return p;
161}
162
163/*
164Bool_t MJCut::ReadTasks(const char *fname, MTask* &env1, MTask* &env2) const
165{
166 // const TString fname = Form("%s/calib%08d.root", fPathIn.Data(), fSequence.GetSequence());
167
168 *fLog << inf << "Reading from file: " << fname << endl;
169
170 TFile file(fname, "READ");
171 if (!file.IsOpen())
172 {
173 *fLog << err << dbginf << "ERROR - Could not open file " << fname << endl;
174 return kFALSE;
175 }
176
177 TObject *o = file.Get("EstimateEnergy");
178 if (o && !o->InheritsFrom(MTask::Class()))
179 {
180 *fLog << err << dbginf << "ERROR - EstimateEnergy read from " << fname << " doesn't inherit from MTask!" << endl;
181 return kFALSE;
182 }
183 env1 = o ? (MTask*)o->Clone() : NULL;
184
185 o = file.Get("CalcHadronness");
186 if (o && !o->InheritsFrom(MTask::Class()))
187 {
188 *fLog << err << dbginf << "ERROR - CalcHadronness read from " << fname << " doesn't inherit from MTask!" << endl;
189 return kFALSE;
190 }
191 env2 = o ? (MTask*)o->Clone() : NULL;
192
193 return kTRUE;
194}
195*/
196
197// --------------------------------------------------------------------------
198//
199// Write the tasks in cont to the file corresponding to analysis number num,
200// see GetOutputFile()
201//
202Bool_t MJCut::WriteTasks(UInt_t num, TObjArray &cont) const
203{
204 if (fPathOut.IsNull())
205 {
206 *fLog << inf << "No output path specified via SetPathOut - no output written." << endl;
207 return kTRUE;
208 }
209
210 const TString oname(GetOutputFile(num));
211
212 *fLog << inf << "Writing to file: " << oname << endl;
213
214 TFile file(oname, fOverwrite?"RECREATE":"NEW", "File created by MJCut", 9);
215 if (!file.IsOpen())
216 {
217 *fLog << err << "ERROR - Couldn't open file " << oname << " for writing..." << endl;
218 return kFALSE;
219 }
220
221 return WriteContainer(cont);
222}
223
224// --------------------------------------------------------------------------
225//
226// Write the result plots and other results to the file corresponding to
227// analysis number num, see GetOutputFile()
228//
229Bool_t MJCut::WriteResult(UInt_t num) const
230{
231 TObjArray arr;
232 return WriteContainer(arr, GetOutputFile(num), "UPDATE");
233}
234
235// --------------------------------------------------------------------------
236//
237// MJCut allows to setup several option by a resource file:
238// MJCut.WriteSummary: yes, no
239// MJCut.SummaryFile: filename
240// MJCut.WriteResult: yes, no
241// MJCut.ResultFile: filename
242//
243Bool_t MJCut::CheckEnvLocal()
244{
245 const TString f0(GetEnv("SummaryFile", ""));
246 const TString f1(GetEnv("ResultFile", ""));
247 if (!f0.IsNull())
248 SetNameSummaryFile(f0);
249 if (!f1.IsNull())
250 SetNameResultFile(f1);
251
252 EnableStorageOfSummary(GetEnv("SummaryFile", fStoreSummary));
253 EnableStorageOfResult(GetEnv("ResultFile", fStoreResult));
254 EnableWobbleMode(GetEnv("WobbleMode", fIsWobble));
255
256 return kTRUE;
257}
258
259// --------------------------------------------------------------------------
260//
261// Setup write to write:
262// container tree optional?
263// -------------- ---------- -----------
264// "MHillas" to "Events"
265// "MHillasSrc" to "Events"
266// "MHadronness" to "Events" yes
267// "MEnergyEst" to "Events" yes
268// "DataType" to "Events"
269//
270void MJCut::SetupWriter(MWriteRootFile *write, const char *name) const
271{
272 if (!write)
273 return;
274
275 write->SetName(name);
276 write->AddContainer("MHillas", "Events");
277 write->AddContainer("MHillasSrc", "Events");
278 write->AddContainer("MHillasSrcAnti", "Events", kFALSE);
279 write->AddContainer("MNewImagePar", "Events", kFALSE);
280 write->AddContainer("MNewImagePar2", "Events", kFALSE);
281 write->AddContainer("MHadronness", "Events", kFALSE);
282 write->AddContainer("MEnergyEst", "Events", kFALSE);
283 write->AddContainer("DataType", "Events");
284
285 // Should not be the default: Either as option, or as
286 // setup from resource file
287 // write.AddContainer("MHillasExt", "Events");
288 // write.AddContainer("MImagePar", "Events");
289 // write.AddContainer("MNewImagePar", "Events");
290}
291
292Bool_t MJCut::ProcessFile(const MDataSet &set)
293{
294 if (!set.IsValid())
295 {
296 *fLog << err << "ERROR - DataSet invalid!" << endl;
297 return kFALSE;
298 }
299
300 CheckEnv();
301
302 // --------------------------------------------------------------------------------
303
304 *fLog << inf;
305 fLog->Separator(GetDescriptor());
306 *fLog << "Perform cuts for data set " << set.GetName() << endl;
307 *fLog << endl;
308
309 // --------------------------------------------------------------------------------
310
311 // Setup Parlist
312 MParList plist;
313 plist.AddToList(this); // take care of fDisplay!
314
315 MParameterI par("DataType");
316 plist.AddToList(&par);
317
318 // Setup Tasklist
319 MTaskList tlist;
320 plist.AddToList(&tlist);
321
322 // La Palma Magic1
323 MObservatory obs;
324 plist.AddToList(&obs);
325
326 // Possible source position (eg. Wobble Mode)
327 MPointingPos source("MSourcePos");
328 if (set.GetSourcePos(source))
329 {
330 plist.AddToList(&source);
331 *fLog << inf << "Using Source Position: " << source.GetTitle() << endl;
332 }
333 else
334 *fLog << inf << "No source position applied..." << endl;
335
336 // Initialize default binnings
337 MBinning bins1(18, 0, 90, "BinningAlpha", "lin");
338 MBinning bins2(15, 10, 1e6 , "BinningEnergyEst", "log");
339 MBinning bins3(50, 0, 60, "BinningTheta", "cos");
340 MBinning bins4("BinningFalseSource");
341 //MBinning bins5("BinningWidth");
342 //MBinning bins6("BinningLength");
343 //MBinning bins7("BinningDist");
344 plist.AddToList(&bins1);
345 plist.AddToList(&bins2);
346 plist.AddToList(&bins3);
347 plist.AddToList(&bins4);
348 //plist.AddToList(&bins5);
349 //plist.AddToList(&bins6);
350 //plist.AddToList(&bins7);
351
352 // --------------------------------------------------------------------------------
353
354 // Setup fitter and histograms
355 MAlphaFitter fit;
356 plist.AddToList(&fit);
357 if (fIsWobble)
358 fit.SetScaleMode(MAlphaFitter::kNone);
359
360 MFillH falpha("MHAlphaOff [MHAlpha]", "MHillasSrc", "FillAlpha");
361 MFillH ffs("MHFalseSourceOff [MHFalseSource]", "MHillas", "FillFS");
362
363 // FIXME: If fPathIn read cuts and energy estimator from file!
364 MContinue cont0("", "Cut0");
365 MContinue cont1("", "Cut1");
366 MContinue cont2("", "Cut2");
367 MContinue cont3("", "Cut3");
368 cont0.SetAllowEmpty();
369 cont1.SetAllowEmpty();
370 cont2.SetAllowEmpty();
371 cont3.SetAllowEmpty();
372
373 // ------------- Loop Off Data --------------------
374 MReadReports readoff;
375 readoff.AddTree("Events", "MTime.", kTRUE);
376 readoff.AddTree("Drive");
377 readoff.AddTree("EffectiveOnTime");
378 if (fIsWobble)
379 set.AddFilesOn(readoff);
380 else
381 set.AddFilesOff(readoff);
382
383 const TString path(Form("%s/", fPathOut.Data()));
384 TString fname0(path);
385 TString fname1(path);
386 fname0 += fNameSummary.IsNull() ? (TString) Form("ganymed%08d-summary.root", set.GetNumAnalysis()) : fNameSummary;
387 fname1 += fNameResult.IsNull() ? (TString) Form("ganymed%08d-result.root", set.GetNumAnalysis()) : fNameResult;
388
389 MWriteRootFile *write0 = CanStoreSummary() ? new MWriteRootFile(fPathOut.IsNull()?0:fname0.Data(), fOverwrite?"RECREATE":"NEW") : 0;
390 MWriteRootFile *write1 = CanStoreResult() ? new MWriteRootFile(fPathOut.IsNull()?0:fname1.Data(), fOverwrite?"RECREATE":"NEW") : 0;
391 SetupWriter(write0, "WriteAfterCut0");
392 SetupWriter(write1, "WriteAfterCut3");
393
394
395 MEnergyEstimate est;
396
397 MTaskEnv taskenv1("EstimateEnergy");
398 taskenv1.SetDefault(fEstimateEnergy ? fEstimateEnergy : &est);
399
400 MTaskEnv taskenv2("CalcHadronness");
401 taskenv2.SetDefault(fCalcHadronness);
402
403 MFillH fill1a("MHHillasOffPre [MHHillas]", "MHillas", "FillHillasPre");
404 MFillH fill2a("MHHillasOffPost [MHHillas]", "MHillas", "FillHillasPost");
405 MFillH fill3a("MHVsSizeOffPost [MHVsSize]", "MHillasSrc", "FillVsSizePost");
406 MFillH fill4a("MHHilExtOffPost [MHHillasExt]", "MHillasSrc", "FillHilExtPost");
407 MFillH fill5a("MHHilSrcOffPost [MHHillasSrc]", "MHillasSrc", "FillHilSrcPost");
408 MFillH fill6a("MHImgParOffPost [MHImagePar]", "MImagePar", "FillImgParPost");
409 MFillH fill7a("MHNewParOffPost [MHNewImagePar]", "MNewImagePar", "FillNewParPost");
410 fill1a.SetNameTab("PreCut");
411 fill2a.SetNameTab("PostCut");
412 fill3a.SetNameTab("VsSize");
413 fill4a.SetNameTab("HilExt");
414 fill5a.SetNameTab("HilSrc");
415 fill6a.SetNameTab("ImgPar");
416 fill7a.SetNameTab("NewPar");
417
418 MPrint print2("MEffectiveOnTime");
419
420 // How to get source position from off- and on-data?
421 MSrcPosCalc scalc;
422 if (fIsWobble)
423 scalc.SetWobbleMode(); /********************/
424 MHillasCalc hcalc;
425 MHillasCalc hcalc2("MHillasCalcAnti");
426 hcalc.SetFlags(MHillasCalc::kCalcHillasSrc);
427 hcalc2.SetFlags(MHillasCalc::kCalcHillasSrc);
428 hcalc2.SetNameHillasSrc("MHillasSrcAnti");
429 hcalc2.SetNameSrcPosCam("MSrcPosAnti");
430
431 MTaskList tlist2;
432 tlist2.AddToList(&scalc);
433 tlist2.AddToList(&hcalc);
434 if (fIsWobble)
435 tlist2.AddToList(&hcalc2);
436 tlist2.AddToList(&taskenv1);
437 tlist2.AddToList(&taskenv2);
438 tlist2.AddToList(&cont0);
439 if (write0)
440 tlist2.AddToList(write0);
441 if (!fWriteOnly)
442 tlist2.AddToList(&fill1a);
443 tlist2.AddToList(&cont1);
444 if (!fWriteOnly && !fIsWobble)
445 tlist2.AddToList(&ffs);
446 tlist2.AddToList(&cont2);
447 if (!fWriteOnly)
448 {
449 tlist2.AddToList(&fill2a);
450 if (fFullDisplay)
451 {
452 tlist2.AddToList(&fill3a);
453 tlist2.AddToList(&fill4a);
454 tlist2.AddToList(&fill5a);
455 tlist2.AddToList(&fill6a);
456 tlist2.AddToList(&fill7a);
457 }
458 }
459 if (!fWriteOnly)
460 tlist2.AddToList(&falpha);
461 tlist2.AddToList(&cont3);
462 if (write1)
463 tlist2.AddToList(write1);
464
465 tlist.AddToList(&readoff);
466 tlist.AddToList(&print2, "EffectiveOnTime");
467 tlist.AddToList(&tlist2, "Events");
468
469 par.SetVal(0);
470
471 // Create and setup the eventloop
472 MEvtLoop evtloop(fName);
473 evtloop.SetParList(&plist);
474 evtloop.SetDisplay(fDisplay);
475 evtloop.SetLogStream(fLog);
476 if (!SetupEnv(evtloop))
477 return kFALSE;
478
479 if (set.HasOffSequences() || fIsWobble)
480 {
481 // Execute first analysis
482 if (!evtloop.Eventloop(fMaxEvents))
483 {
484 *fLog << err << GetDescriptor() << ": Processing of off-sequences failed." << endl;
485 return kFALSE;
486 }
487
488 tlist.PrintStatistics();
489
490 if (!evtloop.GetDisplay())
491 {
492 *fLog << err << GetDescriptor() << ": Execution stopped by user." << endl;
493 return kFALSE;
494 }
495
496 //plist.FindObject("MTimeEffectiveOnTime")->Clear();
497 }
498
499 // ------------- Loop On Data --------------------
500 MReadReports readon;
501 readon.AddTree("Events", "MTime.", kTRUE);
502 readon.AddTree("Drive");
503 readon.AddTree("EffectiveOnTime");
504 set.AddFilesOn(readon);
505
506 if (fIsWobble)
507 scalc.SetWobbleMode(kFALSE); /********************/
508
509 MFillH fill1b("MHHillasOnPre [MHHillas]", "MHillas", "FillHillasPre");
510 MFillH fill2b("MHHillasOnPost [MHHillas]", "MHillas", "FillHillasPost");
511 MFillH fill3b("MHVsSizeOnPost [MHVsSize]", "MHillasSrc", "FillVsSizePost");
512 MFillH fill4b("MHHilExtOnPost [MHHillasExt]", "MHillasSrc", "FillHilExtPost");
513 MFillH fill5b("MHHilSrcOnPost [MHHillasSrc]", "MHillasSrc", "FillHilSrcPost");
514 MFillH fill6b("MHImgParOnPost [MHImagePar]", "MImagePar", "FillImgParPost");
515 MFillH fill7b("MHNewParOnPost [MHNewImagePar]", "MNewImagePar", "FillNewParPost");
516 fill1b.SetNameTab("PreCut");
517 fill2b.SetNameTab("PostCut");
518 fill3b.SetNameTab("VsSize");
519 fill4b.SetNameTab("HilExt");
520 fill5b.SetNameTab("HilSrc");
521 fill6b.SetNameTab("ImgPar");
522 fill7b.SetNameTab("NewPar");
523 fill1b.SetDrawOption(set.HasOffSequences()||fIsWobble?"same":"");
524 fill2b.SetDrawOption(set.HasOffSequences()||fIsWobble?"same":"");
525 fill3b.SetDrawOption(set.HasOffSequences()||fIsWobble?"same":"");
526 fill4b.SetDrawOption(set.HasOffSequences()||fIsWobble?"same":"");
527 fill5b.SetDrawOption(set.HasOffSequences()||fIsWobble?"same":"");
528 fill6b.SetDrawOption(set.HasOffSequences()||fIsWobble?"same":"");
529 fill7b.SetDrawOption(set.HasOffSequences()||fIsWobble?"same":"");
530
531 /*
532 MHVsTime hvs("MEffectiveOnTime.fVal");
533 hvs.SetTitle("Effective On-Time vs. Time;;T_{on}");
534 MFillH fillvs(&hvs, "MTimeEffectiveOnTime", "FillOnTime");
535 fillvs.SetNameTab("OnTime");
536 */
537 MH3 hvs("MPointingPos.fZd");
538 hvs.SetName("Theta");
539 hvs.SetTitle("Effective On-Time vs. Zenith Angle;\\Theta [\\circ];T_{on} [s]");
540 MFillH fillvs(&hvs, "", "FillOnTime");
541 fillvs.SetWeight("MEffectiveOnTime");
542 fillvs.SetNameTab("OnTime");
543
544/*
545 MWeight weight;
546 weight.SetWeight(-1);
547 fill2a.SetWeight(&weight);
548 fill3a.SetWeight(&weight);
549 */
550 MFillH falpha2("MHAlpha", "MHillasSrc", "FillAlpha");
551 MFillH ffs2("MHFalseSource", "MHillas", "FillFS");
552
553 tlist.Replace(&readon);
554 if (!fWriteOnly)
555 {
556 tlist2.Replace(&fill1b);
557/* if (fIsWobble)
558 {
559 tlist2.AddToListAfter(&fill2b, &fill2a);
560 tlist2.AddToListAfter(&fill3b, &fill3a);
561 }
562 else
563 */
564 tlist2.Replace(&fill2b);
565 if (fFullDisplay)
566 {
567 tlist2.Replace(&fill3b);
568 tlist2.Replace(&fill4b);
569 tlist2.Replace(&fill5b);
570 tlist2.Replace(&fill6b);
571 tlist2.Replace(&fill7b);
572 }
573 tlist2.Replace(&falpha2);
574 if (!fIsWobble)
575 tlist2.Replace(&ffs2);
576
577 tlist.AddToList(&fillvs, "EffectiveOnTime");
578 }
579
580 par.SetVal(1);
581
582 TObjArray cont;
583 cont.Add(&cont0);
584 cont.Add(&cont1);
585 cont.Add(&cont2);
586 cont.Add(&cont3);
587 if (taskenv1.GetTask())
588 cont.Add(taskenv1.GetTask());
589 if (taskenv2.GetTask())
590 cont.Add(taskenv2.GetTask());
591
592 if (!WriteTasks(set.GetNumAnalysis(), cont))
593 return kFALSE;
594
595 // Execute first analysis
596 if (!evtloop.Eventloop(fMaxEvents))
597 {
598 *fLog << err << GetDescriptor() << ": Processing of on-sequences failed." << endl;
599 return kFALSE;
600 }
601
602 if (write0)
603 delete write0;
604 if (write1)
605 delete write1;
606
607 tlist.PrintStatistics();
608
609 // FIXME: Perform fit and plot energy dependant alpha plots
610 // and fit result to new tabs!
611 if (!WriteResult(set.GetNumAnalysis()))
612 return kFALSE;
613
614 *fLog << all << GetDescriptor() << ": Done." << endl;
615 *fLog << endl << endl;
616
617 return kTRUE;
618}
Note: See TracBrowser for help on using the repository browser.