source: tags/Mars-V0.9.1/mjobs/MJCut.cc

Last change on this file was 6941, checked in by tbretz, 19 years ago
*** empty log message ***
File size: 19.8 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), /*fSubstraction(kFALSE),*/
77 fEstimateEnergy(0), 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 EnableFullDisplay(GetEnv("FullDisplay", fFullDisplay));
256 //EnableSubstraction(GetEnv("HistogramSubstraction", fSubstraction));
257
258 return kTRUE;
259}
260
261// --------------------------------------------------------------------------
262//
263// Setup write to write:
264// container tree optional?
265// -------------- ---------- -----------
266// "MHillas" to "Events"
267// "MHillasSrc" to "Events"
268// "MHadronness" to "Events" yes
269// "MEnergyEst" to "Events" yes
270// "DataType" to "Events"
271//
272void MJCut::SetupWriter(MWriteRootFile *write, const char *name) const
273{
274 if (!write)
275 return;
276
277 write->SetName(name);
278 write->AddContainer("MHillas", "Events");
279 write->AddContainer("MHillasSrc", "Events");
280 write->AddContainer("MPointingPos", "Events");
281 write->AddContainer("MHillasSrcAnti", "Events", kFALSE);
282 write->AddContainer("MNewImagePar", "Events", kFALSE);
283 write->AddContainer("MNewImagePar2", "Events", kFALSE);
284 write->AddContainer("MHadronness", "Events", kFALSE);
285 write->AddContainer("MEnergyEst", "Events", kFALSE);
286 write->AddContainer("DataType", "Events");
287
288 // Should not be the default: Either as option, or as
289 // setup from resource file
290 // write.AddContainer("MHillasExt", "Events");
291 // write.AddContainer("MImagePar", "Events");
292 // write.AddContainer("MNewImagePar", "Events");
293}
294
295Bool_t MJCut::ProcessFile(const MDataSet &set)
296{
297 if (!set.IsValid())
298 {
299 *fLog << err << "ERROR - DataSet invalid!" << endl;
300 return kFALSE;
301 }
302
303 CheckEnv();
304
305 // --------------------------------------------------------------------------------
306
307 *fLog << inf;
308 fLog->Separator(GetDescriptor());
309 *fLog << "Perform cuts for data set " << set.GetName() << endl;
310 *fLog << endl;
311
312 // --------------------------------------------------------------------------------
313
314 // Setup Parlist
315 MParList plist;
316 plist.AddToList(this); // take care of fDisplay!
317
318 MParameterI par("DataType");
319 plist.AddToList(&par);
320
321 // Setup Tasklist
322 MTaskList tlist;
323 plist.AddToList(&tlist);
324
325 // La Palma Magic1
326 MObservatory obs;
327 plist.AddToList(&obs);
328
329 // Possible source position (eg. Wobble Mode)
330 MPointingPos source("MSourcePos");
331 if (set.GetSourcePos(source))
332 {
333 plist.AddToList(&source);
334 *fLog << inf << "Using Source Position: " << source.GetTitle() << endl;
335 }
336 else
337 *fLog << inf << "No source position applied..." << endl;
338
339 // Initialize default binnings
340 MBinning bins1(18, 0, 90, "BinningAlpha", "lin");
341 MBinning bins2(15, 10, 1e6 , "BinningEnergyEst", "log");
342 MBinning bins3(50, 0, 60, "BinningTheta", "cos");
343 MBinning bins4("BinningFalseSource");
344 //MBinning bins5("BinningWidth");
345 //MBinning bins6("BinningLength");
346 //MBinning bins7("BinningDist");
347 plist.AddToList(&bins1);
348 plist.AddToList(&bins2);
349 plist.AddToList(&bins3);
350 plist.AddToList(&bins4);
351 //plist.AddToList(&bins5);
352 //plist.AddToList(&bins6);
353 //plist.AddToList(&bins7);
354
355 // --------------------------------------------------------------------------------
356
357 // Setup fitter and histograms
358 MAlphaFitter fit;
359 plist.AddToList(&fit);
360 if (fIsWobble)
361 fit.SetScaleMode(MAlphaFitter::kNone);
362
363 MFillH falpha("MHAlphaOff [MHAlpha]", "MHillasSrc", "FillAlpha");
364 MFillH ffs("MHFalseSourceOff [MHFalseSource]", "MHillas", "FillFS");
365
366 // FIXME: If fPathIn read cuts and energy estimator from file!
367 MContinue cont0("", "Cut0");
368 MContinue cont1("", "Cut1");
369 MContinue cont2("", "Cut2");
370 MContinue cont3("", "Cut3");
371 cont0.SetAllowEmpty();
372 cont1.SetAllowEmpty();
373 cont2.SetAllowEmpty();
374 cont3.SetAllowEmpty();
375
376 // ------------- Loop Off Data --------------------
377 MReadReports readoff;
378 readoff.AddTree("Events", "MTime.", kTRUE);
379 readoff.AddTree("Drive");
380 readoff.AddTree("EffectiveOnTime");
381 if (fIsWobble)
382 set.AddFilesOn(readoff);
383 else
384 set.AddFilesOff(readoff);
385
386 const TString path(Form("%s/", fPathOut.Data()));
387 TString fname0(path);
388 TString fname1(path);
389 fname0 += fNameSummary.IsNull() ? (TString) Form("ganymed%08d-summary.root", set.GetNumAnalysis()) : fNameSummary;
390 fname1 += fNameResult.IsNull() ? (TString) Form("ganymed%08d-result.root", set.GetNumAnalysis()) : fNameResult;
391
392 MWriteRootFile *write0 = CanStoreSummary() ? new MWriteRootFile(fPathOut.IsNull()?0:fname0.Data(), fOverwrite?"RECREATE":"NEW") : 0;
393 MWriteRootFile *write1 = CanStoreResult() ? new MWriteRootFile(fPathOut.IsNull()?0:fname1.Data(), fOverwrite?"RECREATE":"NEW") : 0;
394 SetupWriter(write0, "WriteAfterCut0");
395 SetupWriter(write1, "WriteAfterCut3");
396
397
398 MEnergyEstimate est;
399
400 MTaskEnv taskenv1("EstimateEnergy");
401 taskenv1.SetDefault(fEstimateEnergy ? fEstimateEnergy : &est);
402
403 MTaskEnv taskenv2("CalcHadronness");
404 taskenv2.SetDefault(fCalcHadronness);
405
406 MFillH fill1a("MHHillasOffPre [MHHillas]", "MHillas", "FillHillasPre");
407 MFillH fill2a("MHHillasOffPost [MHHillas]", "MHillas", "FillHillasPost");
408 MFillH fill3a("MHVsSizeOffPost [MHVsSize]", "MHillasSrc", "FillVsSizePost");
409 MFillH fill4a("MHHilExtOffPost [MHHillasExt]", "MHillasSrc", "FillHilExtPost");
410 MFillH fill5a("MHHilSrcOffPost [MHHillasSrc]", "MHillasSrc", "FillHilSrcPost");
411 MFillH fill6a("MHImgParOffPost [MHImagePar]", "MImagePar", "FillImgParPost");
412 MFillH fill7a("MHNewParOffPost [MHNewImagePar]", "MNewImagePar", "FillNewParPost");
413 fill1a.SetNameTab("PreCut");
414 fill2a.SetNameTab("PostCut");
415 fill3a.SetNameTab("VsSize");
416 fill4a.SetNameTab("HilExt");
417 fill5a.SetNameTab("HilSrc");
418 fill6a.SetNameTab("ImgPar");
419 fill7a.SetNameTab("NewPar");
420
421 MPrint print2("MEffectiveOnTime");
422
423 // How to get source position from off- and on-data?
424 MSrcPosCalc scalc;
425 if (fIsWobble)
426 scalc.SetWobbleMode(); /********************/
427 MHillasCalc hcalc;
428 MHillasCalc hcalc2("MHillasCalcAnti");
429 hcalc.SetFlags(MHillasCalc::kCalcHillasSrc);
430 hcalc2.SetFlags(MHillasCalc::kCalcHillasSrc);
431 hcalc2.SetNameHillasSrc("MHillasSrcAnti");
432 hcalc2.SetNameSrcPosCam("MSrcPosAnti");
433
434 MTaskList tlist2;
435 tlist2.AddToList(&scalc);
436 tlist2.AddToList(&hcalc);
437 if (fIsWobble)
438 tlist2.AddToList(&hcalc2);
439 tlist2.AddToList(&taskenv1);
440 tlist2.AddToList(&taskenv2);
441 tlist2.AddToList(&cont0);
442 if (write0)
443 tlist2.AddToList(write0);
444 if (!fWriteOnly)
445 tlist2.AddToList(&fill1a);
446 tlist2.AddToList(&cont1);
447 if (!fWriteOnly && !fIsWobble)
448 tlist2.AddToList(&ffs);
449 tlist2.AddToList(&cont2);
450 if (!fWriteOnly)
451 {
452 tlist2.AddToList(&fill2a);
453 if (fFullDisplay)
454 {
455 tlist2.AddToList(&fill3a);
456 tlist2.AddToList(&fill4a);
457 tlist2.AddToList(&fill5a);
458 tlist2.AddToList(&fill6a);
459 tlist2.AddToList(&fill7a);
460 }
461 }
462 if (!fWriteOnly)
463 tlist2.AddToList(&falpha);
464 tlist2.AddToList(&cont3);
465 if (write1)
466 tlist2.AddToList(write1);
467
468 tlist.AddToList(&readoff);
469 tlist.AddToList(&print2, "EffectiveOnTime");
470 tlist.AddToList(&tlist2, "Events");
471
472 par.SetVal(0);
473
474 // Create and setup the eventloop
475 MEvtLoop evtloop(fName);
476 evtloop.SetParList(&plist);
477 evtloop.SetDisplay(fDisplay);
478 evtloop.SetLogStream(fLog);
479 if (!SetupEnv(evtloop))
480 return kFALSE;
481
482 if (set.HasOffSequences() || fIsWobble)
483 {
484 // Execute first analysis
485 if (!evtloop.Eventloop(fMaxEvents))
486 {
487 *fLog << err << GetDescriptor() << ": Processing of off-sequences failed." << endl;
488 return kFALSE;
489 }
490
491 tlist.PrintStatistics();
492
493 if (!evtloop.GetDisplay())
494 {
495 *fLog << err << GetDescriptor() << ": Execution stopped by user." << endl;
496 return kFALSE;
497 }
498
499 //plist.FindObject("MTimeEffectiveOnTime")->Clear();
500 }
501
502 // ------------- Loop On Data --------------------
503 MReadReports readon;
504 readon.AddTree("Events", "MTime.", kTRUE);
505 readon.AddTree("Drive");
506 readon.AddTree("EffectiveOnTime");
507 set.AddFilesOn(readon);
508
509 if (fIsWobble)
510 scalc.SetWobbleMode(kFALSE); /********************/
511
512 MFillH fill1b("MHHillasOnPre [MHHillas]", "MHillas", "FillHillasPre");
513 MFillH fill2b("MHHillasOnPost [MHHillas]", "MHillas", "FillHillasPost");
514 MFillH fill3b("MHVsSizeOnPost [MHVsSize]", "MHillasSrc", "FillVsSizePost");
515 MFillH fill4b("MHHilExtOnPost [MHHillasExt]", "MHillasSrc", "FillHilExtPost");
516 MFillH fill5b("MHHilSrcOnPost [MHHillasSrc]", "MHillasSrc", "FillHilSrcPost");
517 MFillH fill6b("MHImgParOnPost [MHImagePar]", "MImagePar", "FillImgParPost");
518 MFillH fill7b("MHNewParOnPost [MHNewImagePar]", "MNewImagePar", "FillNewParPost");
519 fill1b.SetNameTab("PreCut");
520 fill2b.SetNameTab("PostCut");
521 fill3b.SetNameTab("VsSize");
522 fill4b.SetNameTab("HilExt");
523 fill5b.SetNameTab("HilSrc");
524 fill6b.SetNameTab("ImgPar");
525 fill7b.SetNameTab("NewPar");
526 fill1b.SetDrawOption(set.HasOffSequences()||fIsWobble?"same":"");
527 fill2b.SetDrawOption(set.HasOffSequences()||fIsWobble?"same":"");
528 fill3b.SetDrawOption(set.HasOffSequences()||fIsWobble?"same":"");
529 fill4b.SetDrawOption(set.HasOffSequences()||fIsWobble?"same":"");
530 fill5b.SetDrawOption(set.HasOffSequences()||fIsWobble?"same":"");
531 fill6b.SetDrawOption(set.HasOffSequences()||fIsWobble?"same":"");
532 fill7b.SetDrawOption(set.HasOffSequences()||fIsWobble?"same":"");
533
534 /*
535 MHVsTime hvs("MEffectiveOnTime.fVal");
536 hvs.SetTitle("Effective On-Time vs. Time;;T_{on}");
537 MFillH fillvs(&hvs, "MTimeEffectiveOnTime", "FillOnTime");
538 fillvs.SetNameTab("OnTime");
539 */
540 MH3 hvs("MPointingPos.fZd");
541 hvs.SetName("Theta");
542 hvs.SetTitle("Effective On-Time vs. Zenith Angle;\\Theta [\\circ];T_{on} [s]");
543 MFillH fillvs(&hvs, "", "FillOnTime");
544 fillvs.SetWeight("MEffectiveOnTime");
545 fillvs.SetNameTab("OnTime");
546
547 /*
548 MParameterD weight;
549 weight.SetVal(-1);
550 fill2a.SetWeight(&weight);
551 fill3a.SetWeight(&weight);
552 fill4a.SetWeight(&weight);
553 fill5a.SetWeight(&weight);
554 fill6a.SetWeight(&weight);
555 fill7a.SetWeight(&weight);
556 if (fSubstraction)
557 {
558 fill2a.SetNameTab("PostCut-");
559 fill3a.SetNameTab("VsSize-");
560 fill4a.SetNameTab("HilExt-");
561 fill5a.SetNameTab("HilSrc-");
562 fill6a.SetNameTab("ImgPar-");
563 fill7a.SetNameTab("NewPar-");
564 }
565 */
566 MFillH falpha2("MHAlpha", "MHillasSrc", "FillAlpha");
567 MFillH ffs2("MHFalseSource", "MHillas", "FillFS");
568
569 tlist.Replace(&readon);
570 if (!fWriteOnly)
571 {
572 tlist2.Replace(&fill1b);
573/* if (fIsWobble)
574 {
575 tlist2.AddToListAfter(&fill2b, &fill2a);
576 tlist2.AddToListAfter(&fill3b, &fill3a);
577 }
578 else
579 */
580 tlist2.Replace(&fill2b);
581 if (fFullDisplay)
582 {
583 tlist2.Replace(&fill3b);
584 tlist2.Replace(&fill4b);
585 tlist2.Replace(&fill5b);
586 tlist2.Replace(&fill6b);
587 tlist2.Replace(&fill7b);
588 }
589 tlist2.Replace(&falpha2);
590 if (!fIsWobble)
591 tlist2.Replace(&ffs2);
592
593 tlist.AddToList(&fillvs, "EffectiveOnTime");
594 }
595
596 par.SetVal(1);
597
598 TObjArray cont;
599 cont.Add(&cont0);
600 cont.Add(&cont1);
601 cont.Add(&cont2);
602 cont.Add(&cont3);
603 if (taskenv1.GetTask())
604 cont.Add(taskenv1.GetTask());
605 if (taskenv2.GetTask())
606 cont.Add(taskenv2.GetTask());
607
608 if (!WriteTasks(set.GetNumAnalysis(), cont))
609 return kFALSE;
610
611 // Execute first analysis
612 if (!evtloop.Eventloop(fMaxEvents))
613 {
614 *fLog << err << GetDescriptor() << ": Processing of on-sequences failed." << endl;
615 return kFALSE;
616 }
617
618 if (write0)
619 delete write0;
620 if (write1)
621 delete write1;
622
623 tlist.PrintStatistics();
624
625 // FIXME: Perform fit and plot energy dependant alpha plots
626 // and fit result to new tabs!
627 if (!WriteResult(set.GetNumAnalysis()))
628 return kFALSE;
629
630 *fLog << all << GetDescriptor() << ": Done." << endl;
631 *fLog << endl << endl;
632
633 return kTRUE;
634}
Note: See TracBrowser for help on using the repository browser.