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

Last change on this file since 6437 was 6282, checked in by tbretz, 20 years ago
*** empty log message ***
File size: 14.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 "MBinning.h"
58#include "MSequences.h"
59#include "MParameters.h"
60#include "MObservatory.h"
61
62ClassImp(MJCut);
63
64using namespace std;
65
66// --------------------------------------------------------------------------
67//
68// Default constructor.
69//
70MJCut::MJCut(const char *name, const char *title)
71 : fStoreSummary(kFALSE), fStoreResult(kFALSE), fWriteOnly(kFALSE),
72 fEstimateEnergy(0), fCalcHadronness(0)
73{
74 fName = name ? name : "MJCut";
75 fTitle = title ? title : "Standard program to perform g/h-seperation cuts";
76}
77
78MJCut::~MJCut()
79{
80 if (fEstimateEnergy)
81 delete fEstimateEnergy;
82 if (fCalcHadronness)
83 delete fCalcHadronness;
84}
85
86// --------------------------------------------------------------------------
87//
88// Set the name of the summary file (events after cut0)
89// If you give a name the storage of this file is enabled implicitly.
90// If you give no filename the storage is neither enabled nor disabled,
91// but the storage file name is reset.
92// If no filename is set the default filename is used.
93// You can explicitly enable or disable the storage using EnableStoreOf*()
94// The default argument is no filename.
95//
96void MJCut::SetNameSummaryFile(const char *name)
97{
98 fNameSummary=name;
99 if (!fNameSummary.IsNull())
100 EnableStorageOfSummary();
101}
102
103// --------------------------------------------------------------------------
104//
105// Set the name of the summary file (events after cut3)
106// If you give a name the storage of this file is enabled implicitly.
107// If you give no filename the storage is neither enabled nor disabled,
108// but the storage file name is reset.
109// If no filename is set the default filename is used.
110// You can explicitly enable or disable the storage using EnableStoreOf*()
111// The default argument is no filename.
112//
113void MJCut::SetNameResultFile(const char *name)
114{
115 fNameResult=name;
116 if (!fNameResult.IsNull())
117 EnableStorageOfResult();
118}
119
120// --------------------------------------------------------------------------
121//
122// Setup a task estimating the energy. The given task is cloned.
123//
124void MJCut::SetEnergyEstimator(const MTask *task)
125{
126 if (fEstimateEnergy)
127 delete fEstimateEnergy;
128 fEstimateEnergy = task ? (MTask*)task->Clone() : 0;
129}
130
131// --------------------------------------------------------------------------
132//
133// Setup a task calculating the hadronness. The given task is cloned.
134//
135void MJCut::SetHadronnessCalculator(const MTask *task)
136{
137 if (fCalcHadronness)
138 delete fCalcHadronness;
139 fCalcHadronness = task ? (MTask*)task->Clone() : 0;
140}
141
142// --------------------------------------------------------------------------
143//
144// return fOutputPath+"/ganymed%08d.root", num
145//
146TString MJCut::GetOutputFile(UInt_t num) const
147{
148 TString p(fPathOut);
149 p += "/";
150 p += fNameOutput.IsNull() ? Form("ganymed%08d.root", num) : fNameOutput.Data();
151 return p;
152}
153
154/*
155Bool_t MJCut::ReadTasks(const char *fname, MTask* &env1, MTask* &env2) const
156{
157 // const TString fname = Form("%s/calib%08d.root", fPathIn.Data(), fSequence.GetSequence());
158
159 *fLog << inf << "Reading from file: " << fname << endl;
160
161 TFile file(fname, "READ");
162 if (!file.IsOpen())
163 {
164 *fLog << err << dbginf << "ERROR - Could not open file " << fname << endl;
165 return kFALSE;
166 }
167
168 TObject *o = file.Get("EstimateEnergy");
169 if (o && !o->InheritsFrom(MTask::Class()))
170 {
171 *fLog << err << dbginf << "ERROR - EstimateEnergy read from " << fname << " doesn't inherit from MTask!" << endl;
172 return kFALSE;
173 }
174 env1 = o ? (MTask*)o->Clone() : NULL;
175
176 o = file.Get("CalcHadronness");
177 if (o && !o->InheritsFrom(MTask::Class()))
178 {
179 *fLog << err << dbginf << "ERROR - CalcHadronness read from " << fname << " doesn't inherit from MTask!" << endl;
180 return kFALSE;
181 }
182 env2 = o ? (MTask*)o->Clone() : NULL;
183
184 return kTRUE;
185}
186*/
187
188// --------------------------------------------------------------------------
189//
190// Write the tasks in cont to the file corresponding to analysis number num,
191// see GetOutputFile()
192//
193Bool_t MJCut::WriteTasks(UInt_t num, TObjArray &cont) const
194{
195 if (fPathOut.IsNull())
196 {
197 *fLog << inf << "No output path specified via SetPathOut - no output written." << endl;
198 return kTRUE;
199 }
200
201 const TString oname(GetOutputFile(num));
202
203 *fLog << inf << "Writing to file: " << oname << endl;
204
205 TFile file(oname, fOverwrite?"RECREATE":"NEW", "File created by MJCut", 9);
206 if (!file.IsOpen())
207 {
208 *fLog << err << "ERROR - Couldn't open file " << oname << " for writing..." << endl;
209 return kFALSE;
210 }
211
212 return WriteContainer(cont);
213}
214
215// --------------------------------------------------------------------------
216//
217// Write the result plots and other results to the file corresponding to
218// analysis number num, see GetOutputFile()
219//
220Bool_t MJCut::WriteResult(UInt_t num) const
221{
222 TObjArray arr;
223 return WriteContainer(arr, GetOutputFile(num), "UPDATE");
224}
225
226// --------------------------------------------------------------------------
227//
228// MJCut allows to setup several option by a resource file:
229// MJCut.WriteSummary: yes, no
230// MJCut.SummaryFile: filename
231// MJCut.WriteResult: yes, no
232// MJCut.ResultFile: filename
233//
234Bool_t MJCut::CheckEnvLocal()
235{
236 const TString f0(GetEnv("SummaryFile", ""));
237 const TString f1(GetEnv("ResultFile", ""));
238 if (!f0.IsNull())
239 SetNameSummaryFile(f0);
240 if (!f1.IsNull())
241 SetNameResultFile(f1);
242
243 EnableStorageOfSummary(GetEnv("SummaryFile", fStoreSummary));
244 EnableStorageOfResult(GetEnv("ResultFile", fStoreResult));
245
246 return kTRUE;
247}
248
249// --------------------------------------------------------------------------
250//
251// Setup write to write:
252// container tree optional?
253// -------------- ---------- -----------
254// "MHillas" to "Events"
255// "MHillasSrc" to "Events"
256// "MHadronness" to "Events" yes
257// "DataType" to "Events"
258//
259void MJCut::SetupWriter(MWriteRootFile &write) const
260{
261 write.AddContainer("MHillas", "Events");
262 write.AddContainer("MHillasSrc", "Events");
263 write.AddContainer("MHadronness", "Events", kFALSE);
264 write.AddContainer("DataType", "Events");
265
266 // Should not be the default: Either as option, or as
267 // setup from resource file
268 // write.AddContainer("MHillasExt", "Events");
269 // write.AddContainer("MImagePar", "Events");
270 // write.AddContainer("MNewImagePar", "Events");
271}
272
273Bool_t MJCut::ProcessFile(const MSequences &seq)
274{
275 if (!seq.IsValid())
276 {
277 *fLog << err << "ERROR - Sequences invalid!" << endl;
278 return kFALSE;
279 }
280
281 CheckEnv();
282
283 // --------------------------------------------------------------------------------
284
285 *fLog << inf;
286 fLog->Separator(GetDescriptor());
287 *fLog << "Perform cuts for sequences " << seq.GetName() << endl;
288 *fLog << endl;
289
290 // --------------------------------------------------------------------------------
291
292 // Setup Parlist
293 MParList plist;
294 plist.AddToList(this); // take care of fDisplay!
295
296 MParameterI par("DataType");
297 plist.AddToList(&par);
298
299 // Setup Tasklist
300 MTaskList tlist;
301 plist.AddToList(&tlist);
302
303 // La Palma Magic1
304 MObservatory obs;
305 plist.AddToList(&obs);
306
307 // Initialize binnings
308 MBinning bins1(18, 0, 90, "BinningAlpha", "lin");
309 MBinning bins2(25, 10, 1000000, "BinningEnergyEst", "log");
310 MBinning bins3(50, 0, 60, "BinningTheta", "cos");
311 MBinning bins4("BinningFalseSource");
312 plist.AddToList(&bins1);
313 plist.AddToList(&bins2);
314 plist.AddToList(&bins3);
315 plist.AddToList(&bins4);
316
317 // --------------------------------------------------------------------------------
318
319 // Setup fitter and histograms
320 MAlphaFitter fit;
321 plist.AddToList(&fit);
322
323 MFillH falpha("MHAlphaOff [MHAlpha]", "MHillasSrc", "FillAlpha");
324 MFillH ffs("MHFalseSourceOff [MHFalseSource]", "MHillas", "FillFS");
325
326 // FIXME: If fPathIn read cuts and energy estimator from file!
327 MContinue cont0("", "Cut0");
328 MContinue cont1("", "Cut1");
329 MContinue cont2("", "Cut2");
330 MContinue cont3("", "Cut3");
331 cont0.SetAllowEmpty();
332 cont1.SetAllowEmpty();
333 cont2.SetAllowEmpty();
334 cont3.SetAllowEmpty();
335
336 // ------------- Loop Off Data --------------------
337 MReadReports readoff;
338 readoff.AddTree("Events", "MTime.", kTRUE);
339 readoff.AddTree("Drive");
340 readoff.AddTree("EffectiveOnTime");
341 seq.AddFilesOff(readoff);
342
343 const TString path(Form("%s/", fPathOut.Data()));
344 TString fname0(path);
345 TString fname1(path);
346 fname0 += fNameSummary.IsNull() ? Form("ganymed%08d-summary.root", seq.GetNumAnalysis()) : fNameSummary;
347 fname1 += fNameResult.IsNull() ? Form("ganymed%08d-result.root", seq.GetNumAnalysis()) : fNameResult;
348
349 MWriteRootFile write0(fPathOut.IsNull()?0:fname0.Data(), fOverwrite?"RECREATE":"NEW");
350 MWriteRootFile write1(fPathOut.IsNull()?0:fname1.Data(), fOverwrite?"RECREATE":"NEW");
351 if (CanStoreSummary())
352 SetupWriter(write0);
353 if (CanStoreSummary())
354 SetupWriter(write1);
355
356 MEnergyEstimate est;
357
358 MTaskEnv taskenv1("EstimateEnergy");
359 taskenv1.SetDefault(fEstimateEnergy ? fEstimateEnergy : &est);
360
361 MTaskEnv taskenv2("CalcHadronness");
362 taskenv2.SetDefault(fCalcHadronness);
363
364 MFillH fill1a("MHHillasOffPre [MHHillas]", "MHillas", "FillHillasPre");
365 MFillH fill2a("MHHillasOffPost [MHHillas]", "MHillas", "FillHillasPost");
366 fill1a.SetNameTab("PreCut");
367 fill2a.SetNameTab("PostCut");
368
369 MPrint print2("MEffectiveOnTime");
370
371 // How to get source position from off- and on-data?
372 MSrcPosCalc scalc;
373
374 //MHillasCalc hcalc;
375 //hcalc.SetFlags(MHillasCalc::kCalcHillasSrc);
376
377 MTaskList tlist2;
378 tlist2.AddToList(&taskenv1);
379 tlist2.AddToList(&taskenv2);
380 tlist2.AddToList(&cont0);
381 if (CanStoreSummary())
382 tlist2.AddToList(&write0);
383 if (!fWriteOnly)
384 tlist2.AddToList(&fill1a);
385 tlist2.AddToList(&cont1);
386 if (!fWriteOnly)
387 tlist2.AddToList(&ffs);
388 //tlist2.AddToList(&scalc);
389 //tlist2.AddToList(&hcalc);
390 tlist2.AddToList(&cont2);
391 if (!fWriteOnly)
392 tlist2.AddToList(&fill2a);
393 if (!fWriteOnly)
394 tlist2.AddToList(&falpha);
395 tlist2.AddToList(&cont3);
396 if (CanStoreResult())
397 tlist2.AddToList(&write1);
398
399 tlist.AddToList(&readoff);
400 tlist.AddToList(&print2, "EffectiveOnTime");
401 tlist.AddToList(&tlist2, "Events");
402
403 par.SetVal(0);
404
405 // Create and setup the eventloop
406 MEvtLoop evtloop(fName);
407 evtloop.SetParList(&plist);
408 evtloop.SetDisplay(fDisplay);
409 evtloop.SetLogStream(fLog);
410 if (!SetupEnv(evtloop))
411 return kFALSE;
412
413 if (seq.HasOffSequences())
414 {
415 // Execute first analysis
416 if (!evtloop.Eventloop(fMaxEvents))
417 {
418 *fLog << err << GetDescriptor() << ": Processing of off-sequences failed." << endl;
419 return kFALSE;
420 }
421
422 tlist.PrintStatistics();
423
424 if (!evtloop.GetDisplay())
425 {
426 *fLog << err << GetDescriptor() << ": Execution stopped by used." << endl;
427 return kFALSE;
428 }
429 }
430
431 // ------------- Loop On Data --------------------
432 MReadReports readon;
433 readon.AddTree("Events", "MTime.", kTRUE);
434 readon.AddTree("Drive");
435 readon.AddTree("EffectiveOnTime");
436 seq.AddFilesOn(readon);
437
438 MFillH fill1b("MHHillasOnPre [MHHillas]", "MHillas", "FillHillasPre");
439 MFillH fill2b("MHHillasOnPost [MHHillas]", "MHillas", "FillHillasPost");
440 fill1b.SetNameTab("PreCut");
441 fill2b.SetNameTab("PostCut");
442 fill1b.SetDrawOption("same");
443 fill2b.SetDrawOption("same");
444
445 MFillH falpha2("MHAlpha", "MHillasSrc", "FillAlpha");
446 MFillH ffs2("MHFalseSource", "MHillas", "FillFS");
447
448 tlist.Replace(&readon);
449 if (!fWriteOnly)
450 {
451 tlist2.Replace(&fill1b);
452 tlist2.Replace(&fill2b);
453 tlist2.Replace(&falpha2);
454 tlist2.Replace(&ffs2);
455 }
456
457 par.SetVal(1);
458
459 TObjArray cont;
460 cont.Add(&cont0);
461 cont.Add(&cont1);
462 cont.Add(&cont2);
463 cont.Add(&cont3);
464 if (taskenv1.GetTask())
465 cont.Add(taskenv1.GetTask());
466 if (taskenv2.GetTask())
467 cont.Add(taskenv2.GetTask());
468
469 if (!WriteTasks(seq.GetNumAnalysis(), cont))
470 return kFALSE;
471
472 // Execute first analysis
473 if (!evtloop.Eventloop(fMaxEvents))
474 {
475 *fLog << err << GetDescriptor() << ": Processing of on-sequences failed." << endl;
476 return kFALSE;
477 }
478
479 tlist.PrintStatistics();
480
481 // FIXME: Perform fit and plot energy dependant alpha plots
482 // and fit result to new tabs!
483 if (!WriteResult(seq.GetNumAnalysis()))
484 return kFALSE;
485
486 *fLog << all << GetDescriptor() << ": Done." << endl;
487 *fLog << endl << endl;
488
489 return kTRUE;
490}
Note: See TracBrowser for help on using the repository browser.