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

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