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

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