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

Last change on this file since 7517 was 7517, checked in by tbretz, 19 years ago
*** empty log message ***
File size: 25.5 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// Root
35#include <TEnv.h>
36#include <TFile.h>
37
38// Environment
39#include "MLog.h"
40#include "MLogManip.h"
41
42// Eventloop
43#include "MParList.h"
44#include "MTaskList.h"
45#include "MEvtLoop.h"
46
47// Display
48#include "MStatusDisplay.h"
49
50// Tasks
51#include "MReadReports.h"
52#include "MReadMarsFile.h"
53#include "MPrint.h"
54#include "MContinue.h"
55#include "MTaskEnv.h"
56#include "MPointingDevCalc.h"
57#include "MSrcPosCalc.h"
58#include "MSrcPosCorrect.h"
59#include "MHillasCalc.h"
60#include "MFillH.h"
61#include "MWriteRootFile.h"
62
63// Fit signal environment
64#include "../mhflux/MAlphaFitter.h"
65#include "../mhflux/MHAlpha.h"
66
67// Containers
68#include "MH3.h"
69#include "MBinning.h"
70#include "MDataSet.h"
71#include "MParameters.h"
72#include "MPointingPos.h"
73#include "MObservatory.h"
74
75ClassImp(MJCut);
76
77using namespace std;
78
79// --------------------------------------------------------------------------
80//
81// Default constructor. Set defaults for fStoreSummary, fStoreresult,
82// fWriteOnly, fIsWobble and fFullDisplay to kFALSE and initialize
83// /*fEstimateEnergy and*/ fCalcHadronness with NULL.
84//
85MJCut::MJCut(const char *name, const char *title)
86 : fStoreSummary(kFALSE), fStoreResult(kTRUE), fWriteOnly(kFALSE),
87 fIsWobble(kFALSE), fIsMonteCarlo(kFALSE), fFullDisplay(kTRUE),
88 fNameHist("MHThetaSq"), fCalcHadronness(0), fCalcDisp(0)
89{
90 fName = name ? name : "MJCut";
91 fTitle = title ? title : "Standard program to perform g/h-separation cuts";
92}
93
94// --------------------------------------------------------------------------
95//
96// Destructor. Delete fEstimateEnergy and fCalcHadronness if != NULL
97//
98MJCut::~MJCut()
99{
100 //if (fEstimateEnergy)
101 // delete fEstimateEnergy;
102 if (fCalcHadronness)
103 delete fCalcHadronness;
104 if (fCalcDisp)
105 delete fCalcDisp;
106}
107
108// --------------------------------------------------------------------------
109//
110// Set the name of the summary file (events after cut0)
111// If you give a name the storage of this file is enabled implicitly.
112// If you give no filename the storage is neither enabled nor disabled,
113// but the storage file name is reset.
114// If no filename is set the default filename is used.
115// You can explicitly enable or disable the storage using EnableStoreOf*()
116// The default argument is no filename.
117//
118void MJCut::SetNameSummaryFile(const char *name)
119{
120 fNameSummary=name;
121 if (!fNameSummary.IsNull())
122 EnableStorageOfSummary();
123}
124
125// --------------------------------------------------------------------------
126//
127// Set the name of the summary file (events after cut3)
128// If you give a name the storage of this file is enabled implicitly.
129// If you give no filename the storage is neither enabled nor disabled,
130// but the storage file name is reset.
131// If no filename is set the default filename is used.
132// You can explicitly enable or disable the storage using EnableStoreOf*()
133// The default argument is no filename.
134//
135void MJCut::SetNameResultFile(const char *name)
136{
137 fNameResult=name;
138 if (!fNameResult.IsNull())
139 EnableStorageOfResult();
140}
141
142// --------------------------------------------------------------------------
143//
144// Setup a task estimating the energy. The given task is cloned.
145//
146/*
147void MJCut::SetEnergyEstimator(const MTask *task)
148{
149 if (fEstimateEnergy)
150 delete fEstimateEnergy;
151 fEstimateEnergy = task ? (MTask*)task->Clone() : 0;
152}
153*/
154
155// --------------------------------------------------------------------------
156//
157// Setup a task calculating the hadronness. The given task is cloned.
158//
159void MJCut::SetHadronnessCalculator(const MTask *task)
160{
161 if (fCalcHadronness)
162 delete fCalcHadronness;
163 fCalcHadronness = task ? (MTask*)task->Clone() : 0;
164}
165
166// --------------------------------------------------------------------------
167//
168// Setup a task calculating disp. The given task is cloned.
169//
170void MJCut::SetDispCalculator(const MTask *task)
171{
172 if (fCalcDisp)
173 delete fCalcDisp;
174 fCalcDisp = task ? (MTask*)task->Clone() : 0;
175}
176
177// --------------------------------------------------------------------------
178//
179// return fOutputPath+"/ganymed%08d.root", num
180//
181TString MJCut::GetOutputFile(UInt_t num) const
182{
183 TString p(fPathOut);
184 p += "/";
185 p += fNameOutput.IsNull() ? Form("ganymed%08d.root", num) : fNameOutput.Data();
186 gSystem->ExpandPathName(p);
187 return p;
188}
189
190/*
191Bool_t MJCut::ReadTasks(const char *fname, MTask* &env1, MTask* &env2) const
192{
193 // const TString fname = Form("%s/calib%08d.root", fPathIn.Data(), fSequence.GetSequence());
194
195 *fLog << inf << "Reading from file: " << fname << endl;
196
197 TFile file(fname, "READ");
198 if (!file.IsOpen())
199 {
200 *fLog << err << dbginf << "ERROR - Could not open file " << fname << endl;
201 return kFALSE;
202 }
203
204 TObject *o = file.Get("EstimateEnergy");
205 if (o && !o->InheritsFrom(MTask::Class()))
206 {
207 *fLog << err << dbginf << "ERROR - EstimateEnergy read from " << fname << " doesn't inherit from MTask!" << endl;
208 return kFALSE;
209 }
210 env1 = o ? (MTask*)o->Clone() : FNULL;
211
212 o = file.Get("CalcHadronness");
213 if (o && !o->InheritsFrom(MTask::Class()))
214 {
215 *fLog << err << dbginf << "ERROR - CalcHadronness read from " << fname << " doesn't inherit from MTask!" << endl;
216 return kFALSE;
217 }
218 env2 = o ? (MTask*)o->Clone() : NULL;
219
220 return kTRUE;
221}
222*/
223
224// --------------------------------------------------------------------------
225//
226// Write the tasks in cont to the file corresponding to analysis number num,
227// see GetOutputFile()
228//
229Bool_t MJCut::WriteTasks(UInt_t num, TObjArray &cont) const
230{
231 if (fPathOut.IsNull())
232 {
233 *fLog << inf << "No output path specified via SetPathOut - no output written." << endl;
234 return kTRUE;
235 }
236
237 const TString oname(GetOutputFile(num));
238
239 *fLog << inf << "Writing to file: " << oname << endl;
240
241 TFile *file = 0;
242 if (fNameResult.IsNull() && fStoreResult)
243 {
244 file = (TFile*)gROOT->GetListOfFiles()->FindObject(oname);
245 if (file)
246 file->cd();
247 }
248 else
249 file = TFile::Open(oname, fOverwrite?"RECREATE":"NEW", "File created by MJCut", 9);
250
251 if (!file)
252 {
253 *fLog << err << "ERROR - Couldn't open file " << oname << " for writing..." << endl;
254 return kFALSE;
255 }
256
257 const Bool_t rc = WriteContainer(cont);
258
259 if (!(fNameResult.IsNull() && fStoreResult))
260 delete file;
261
262 return rc;
263}
264
265// --------------------------------------------------------------------------
266//
267// Write the result plots and other results to the file corresponding to
268// analysis number num, see GetOutputFile()
269//
270Bool_t MJCut::WriteResult(const MParList &plist, UInt_t num) const
271{
272 TObjArray arr;
273
274 // Save all MBinnings
275 TIter Next(plist);
276 TObject *o=0;
277 while ((o=Next()))
278 if (o->InheritsFrom(MBinning::Class()))
279 arr.Add(o);
280
281 // Save also the result, not only the setup
282 const MHAlpha *halpha = (MHAlpha*)plist.FindObject("Hist", "MHAlpha");
283 if (halpha)
284 arr.Add((TObject*)(&halpha->GetAlphaFitter()));
285
286 const TString fname(fNameOutput.IsNull() ? Form("ganymed%08d.root", num) : fNameOutput.Data());
287
288 // If requested, write to already open output file
289 if (fNameResult.IsNull() && fStoreResult)
290 {
291 TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject(fname);
292 if (file)
293 {
294 file->cd();
295 return WriteContainer(arr);
296 }
297 }
298
299 return WriteContainer(arr, fname, "UPDATE");
300}
301
302// --------------------------------------------------------------------------
303//
304// MJCut allows to setup several option by a resource file:
305// MJCut.WriteSummary: yes, no
306// MJCut.SummaryFile: filename
307// MJCut.WriteResult: yes, no
308// MJCut.ResultFile: filename
309// MJCut.HistName: MHAlpha
310//
311Bool_t MJCut::CheckEnvLocal()
312{
313 const TString f0(GetEnv("SummaryFile", ""));
314 const TString f1(GetEnv("ResultFile", ""));
315 if (!f0.IsNull())
316 SetNameSummaryFile(f0);
317 if (!f1.IsNull())
318 SetNameResultFile(f1);
319
320 EnableStorageOfSummary(GetEnv("SummaryFile", fStoreSummary));
321 EnableStorageOfResult(GetEnv("ResultFile", fStoreResult));
322 EnableWobbleMode(GetEnv("WobbleMode", fIsWobble));
323 EnableMonteCarloMode(GetEnv("MonteCarlo", fIsMonteCarlo));
324 EnableFullDisplay(GetEnv("FullDisplay", fFullDisplay));
325 //EnableSubstraction(GetEnv("HistogramSubstraction", fSubstraction));
326
327 SetNameHist(GetEnv("NameHist", fNameHist));
328 SetNameHistFS(GetEnv("NameHistFS", fNameHistFS));
329
330 return kTRUE;
331}
332
333// --------------------------------------------------------------------------
334//
335// Setup write to write:
336// container tree optional?
337// -------------- ---------- -----------
338// "MHillas" to "Events"
339// "MHillasSrc" to "Events"
340// "Hadronness" to "Events" yes
341// "MEnergyEst" to "Events" yes
342// "DataType" to "Events"
343//
344void MJCut::SetupWriter(MWriteRootFile *write, const char *name) const
345{
346 if (!write)
347 return;
348
349 write->SetName(name);
350 write->AddContainer("MHillas", "Events");
351 write->AddContainer("MHillasSrc", "Events");
352 write->AddContainer("MHillasExt", "Events");
353 write->AddContainer("MPointingPos", "Events");
354 write->AddContainer("MHillasSrcAnti", "Events", kFALSE);
355 write->AddContainer("MImagePar", "Events", kFALSE);
356 write->AddContainer("MNewImagePar", "Events", kFALSE);
357 write->AddContainer("MNewImagePar2", "Events", kFALSE);
358 write->AddContainer("Hadronness", "Events", kFALSE);
359 write->AddContainer("MSrcPosCam", "Events", kFALSE);
360 write->AddContainer("MSrcPosAnti", "Events", kFALSE);
361 write->AddContainer("ThetaSquared", "Events", kFALSE);
362 write->AddContainer("OpticalAxis", "Events", kFALSE);
363 write->AddContainer("Disp", "Events", kFALSE);
364 write->AddContainer("MTime", "Events", kFALSE);
365 write->AddContainer("MMcEvt", "Events", kFALSE);
366 write->AddContainer("DataType", "Events");
367 // write->AddContainer("MMuonSearchPar", "Events", kFALSE);
368 // write->AddContainer("MMuonCalibPar", "Events", kFALSE);
369}
370
371// --------------------------------------------------------------------------
372//
373// Create a new instance of an object with name name of class
374// type fNameHist in parlist. It must derive from MHAlpha.
375// Call ForceUsingSize for it and return its pointer.
376// If something fails NULL is returned.
377//
378MHAlpha *MJCut::CreateNewHist(MParList &plist, const char *name) const
379{
380 TClass *cls = gROOT->GetClass(fNameHist);
381 if (!cls)
382 {
383 *fLog << err << "Class " << fNameHist << " not found in dictionary... abort." << endl;
384 return NULL;
385 }
386 if (!cls->InheritsFrom(MHAlpha::Class()))
387 {
388 *fLog << err << "Class " << fNameHist << " doesn't inherit from MHAlpha... abort." << endl;
389 return NULL;
390 }
391
392 const TString objname(Form("Hist%s", name));
393 MHAlpha *h = (MHAlpha*)plist.FindCreateObj(fNameHist, objname);
394 if (!h)
395 return NULL;
396
397 h->ForceUsingSize();
398
399 return h;
400}
401
402// --------------------------------------------------------------------------
403//
404// Create a new instance of an object with name name of class
405// type fNameHistFS in parlist. It must derive from MHFalseSource
406// If something fails NULL is returned.
407//
408MH *MJCut::CreateNewHistFS(MParList &plist, const char *name) const
409{
410 const TString cname(fNameHistFS.IsNull()?"MHFalseSource":fNameHistFS);
411
412 TClass *cls = gROOT->GetClass(cname);
413 if (!cls)
414 {
415 *fLog << err << "Class " << cname << " not found in dictionary... abort." << endl;
416 return NULL;
417 }
418 if (!cls->InheritsFrom("MHFalseSource"))
419 {
420 *fLog << err << "Class " << cname << " doesn't inherit from MHFalseSource... abort." << endl;
421 return NULL;
422 }
423
424 const TString objname(Form("FS%s", name));
425 return (MH*)plist.FindCreateObj(cname, objname);
426}
427
428Bool_t MJCut::Process(const MDataSet &set)
429{
430 if (!set.IsValid())
431 {
432 *fLog << err << "ERROR - DataSet invalid!" << endl;
433 return kFALSE;
434 }
435
436 CheckEnv();
437
438 // --------------------------------------------------------------------------------
439
440 *fLog << inf;
441 fLog->Separator(GetDescriptor());
442 *fLog << "Perform cuts for data set " << set.GetName() << endl;
443 *fLog << endl;
444
445 // --------------------------------------------------------------------------------
446
447 // Setup Parlist
448 MParList plist;
449 plist.AddToList(this); // take care of fDisplay!
450
451 MParameterI par("DataType");
452 plist.AddToList(&par);
453
454 // Setup Tasklist
455 MTaskList tlist;
456 plist.AddToList(&tlist);
457
458 // La Palma Magic1
459 MObservatory obs;
460 plist.AddToList(&obs);
461
462 // Possible source position (eg. Wobble Mode)
463 MPointingPos source("MSourcePos");
464 if (set.HasSource())
465 {
466 if (!set.GetSourcePos(source))
467 return kFALSE;
468 plist.AddToList(&source);
469 *fLog << all;
470 source.Print("RaDec");
471 }
472 else
473 *fLog << inf << "No source position applied..." << endl;
474
475 // Initialize default binnings
476 MBinning bins1(18, 0, 90, "BinningAlpha", "lin");
477 MBinning bins2(15, 10, 1e6 , "BinningSize", "log");
478 MBinning bins3(67, -0.005, 0.665, "BinningTheta", "asin");
479 MBinning bins4("BinningFalseSource");
480 MBinning bins5("BinningWidth");
481 MBinning bins6("BinningLength");
482 MBinning bins7("BinningDist");
483 MBinning bins8("BinningMaxDist");
484 MBinning bins9("BinningM3Long");
485 MBinning bins0("BinningConc1");
486 plist.AddToList(&bins1);
487 plist.AddToList(&bins2);
488 plist.AddToList(&bins3);
489 plist.AddToList(&bins4);
490 plist.AddToList(&bins5);
491 plist.AddToList(&bins6);
492 plist.AddToList(&bins7);
493 plist.AddToList(&bins8);
494 plist.AddToList(&bins9);
495 plist.AddToList(&bins0);
496 //plist.AddToList(&binsa);
497
498 // --------------------------------------------------------------------------------
499
500 // Setup fitter and histograms
501 MAlphaFitter fit;
502 plist.AddToList(&fit);
503 if (fIsWobble)
504 fit.SetScaleMode(MAlphaFitter::kNone);
505
506 MHAlpha *halphaoff = CreateNewHist(plist, "Off");
507 MFillH falpha(halphaoff, "", "FillHist");
508 MH *hfsoff = CreateNewHistFS(plist, "Off");
509 MFillH ffs(hfsoff, "MHillas", "FillFS");
510
511 // FIXME: If fPathIn read cuts and energy estimator from file!
512 MContinue cont0("", "Cut0");
513 MContinue cont1("", "Cut1");
514 MContinue cont2("", "Cut2");
515 MContinue cont3("", "Cut3");
516 cont0.SetAllowEmpty();
517 cont1.SetAllowEmpty();
518 cont2.SetAllowEmpty();
519 cont3.SetAllowEmpty();
520
521 // ------------- Loop Off Data --------------------
522 MReadReports readoffdata;
523 readoffdata.AddTree("Events", "MTime.", MReadReports::kMaster);
524 readoffdata.AddTree("Drive", MReadReports::kRequired);
525 readoffdata.AddTree("Starguider", MReadReports::kRequired);
526 readoffdata.AddTree("EffectiveOnTime");
527
528 MReadMarsFile readoffmc("Events");
529 readoffmc.DisableAutoScheme();
530
531 MRead &readoff = fIsMonteCarlo ? (MRead&)readoffmc : (MRead&)readoffdata;
532 if (fIsWobble)
533 set.AddFilesOn(readoff);
534 else
535 set.AddFilesOff(readoff);
536
537 const TString path(Form("%s/", fPathOut.Data()));
538 TString fname0(path);
539 TString fname1(path);
540 fname0 += fNameSummary.IsNull() ? (TString) Form("ganymed%08d-summary.root", set.GetNumAnalysis()) : fNameSummary;
541 fname1 += fNameResult.IsNull() ? (TString) Form("ganymed%08d.root", set.GetNumAnalysis()) : fNameResult;
542
543 MWriteRootFile *write0 = CanStoreSummary() ? new MWriteRootFile(fPathOut.IsNull()?0:fname0.Data(), fOverwrite?"RECREATE":"NEW") : 0;
544 MWriteRootFile *write1 = CanStoreResult() ? new MWriteRootFile(fPathOut.IsNull()?0:fname1.Data(), fOverwrite?"RECREATE":"NEW") : 0;
545 SetupWriter(write0, "WriteAfterCut0");
546 SetupWriter(write1, "WriteAfterCut3");
547
548/*
549 MEnergyEstimate est;
550
551 MTaskEnv taskenv1("EstimateEnergy");
552 taskenv1.SetDefault(fEstimateEnergy ? fEstimateEnergy : &est);
553 */
554 MTaskEnv taskenv2("CalcHadronness");
555 taskenv2.SetDefault(fCalcHadronness);
556
557 MTaskEnv taskenv3("CalcDisp");
558 taskenv3.SetDefault(fCalcDisp);
559
560 MFillH fill1a("MHHillasOffPre [MHHillas]", "MHillas", "FillHillasPre");
561 MFillH fill2a("MHHillasOffPost [MHHillas]", "MHillas", "FillHillasPost");
562 MFillH fill3a("MHVsSizeOffPost [MHVsSize]", "MHillasSrc", "FillVsSizePost");
563 MFillH fill4a("MHHilExtOffPost [MHHillasExt]", "MHillasSrc", "FillHilExtPost");
564 MFillH fill5a("MHHilSrcOffPost [MHHillasSrc]", "MHillasSrc", "FillHilSrcPost");
565 MFillH fill6a("MHImgParOffPost [MHImagePar]", "MImagePar", "FillImgParPost");
566 MFillH fill7a("MHNewParOffPost [MHNewImagePar]", "MNewImagePar", "FillNewParPost");
567 fill1a.SetNameTab("PreCut");
568 fill2a.SetNameTab("PostCut");
569 fill3a.SetNameTab("VsSize");
570 fill4a.SetNameTab("HilExt");
571 fill5a.SetNameTab("HilSrc");
572 fill6a.SetNameTab("ImgPar");
573 fill7a.SetNameTab("NewPar");
574
575 MPrint print2("MEffectiveOnTime");
576 print2.EnableSkip();
577
578 // How to get source position from off- and on-data?
579 MSrcPosCalc scalc;
580 scalc.SetMode(fIsWobble?MSrcPosCalc::kWobble:MSrcPosCalc::kOffData); /********************/
581
582 MSrcPosCorrect scor;
583
584 MHillasCalc hcalc;
585 MHillasCalc hcalc2("MHillasCalcAnti");
586 hcalc.SetFlags(MHillasCalc::kCalcHillasSrc);
587 hcalc2.SetFlags(MHillasCalc::kCalcHillasSrc);
588 hcalc2.SetNameHillasSrc("MHillasSrcAnti");
589 hcalc2.SetNameSrcPosCam("MSrcPosAnti");
590
591 MTaskList tlist2;
592 tlist2.AddToList(&scalc);
593 tlist2.AddToList(&scor);
594 tlist2.AddToList(&hcalc);
595 if (fIsWobble)
596 tlist2.AddToList(&hcalc2);
597 //tlist2.AddToList(&taskenv1);
598 tlist2.AddToList(&taskenv2);
599 tlist2.AddToList(&taskenv3);
600 tlist2.AddToList(&cont0);
601 if (write0)
602 tlist2.AddToList(write0);
603 if (!fWriteOnly)
604 tlist2.AddToList(&fill1a);
605 tlist2.AddToList(&cont1);
606 if (!fWriteOnly && (!fIsWobble || !fNameHistFS.IsNull()))
607 tlist2.AddToList(&ffs);
608 tlist2.AddToList(&cont2);
609 if (!fWriteOnly)
610 {
611 tlist2.AddToList(&fill2a);
612 if (fFullDisplay)
613 {
614 tlist2.AddToList(&fill3a);
615 tlist2.AddToList(&fill4a);
616 tlist2.AddToList(&fill5a);
617 tlist2.AddToList(&fill6a);
618 tlist2.AddToList(&fill7a);
619 }
620 }
621 if (!fWriteOnly)
622 tlist2.AddToList(&falpha);
623 tlist2.AddToList(&cont3);
624 if (write1)
625 tlist2.AddToList(write1);
626
627 MPointingDevCalc devcalc;
628
629 tlist.AddToList(&readoff);
630 if (gLog.GetDebugLevel()>4)
631 tlist.AddToList(&print2, "EffectiveOnTime");
632 tlist.AddToList(&devcalc, "Starguider");
633 tlist.AddToList(&tlist2, "Events");
634
635 par.SetVal(0);
636
637 // Create and setup the eventloop
638 MEvtLoop evtloop(fName);
639 evtloop.SetParList(&plist);
640 evtloop.SetDisplay(fDisplay);
641 evtloop.SetLogStream(fLog);
642 if (!SetupEnv(evtloop))
643 return kFALSE;
644
645 TObjArray cont;
646 cont.Add(&cont0);
647 cont.Add(&cont1);
648 cont.Add(&cont2);
649 cont.Add(&cont3);
650 //if (taskenv1.GetTask())
651 // cont.Add(taskenv1.GetTask());
652 if (taskenv2.GetTask())
653 cont.Add(taskenv2.GetTask());
654 if (taskenv3.GetTask())
655 cont.Add(taskenv3.GetTask());
656
657 if (!WriteTasks(set.GetNumAnalysis(), cont))
658 return kFALSE;
659
660 if (set.HasOffSequences() || fIsWobble)
661 {
662 // Execute first analysis
663 if (!evtloop.Eventloop(fMaxEvents))
664 {
665 *fLog << err << GetDescriptor() << ": Processing of off-sequences failed." << endl;
666 return kFALSE;
667 }
668
669 if (!evtloop.GetDisplay())
670 {
671 *fLog << err << GetDescriptor() << ": Execution stopped by user." << endl;
672 return kFALSE;
673 }
674
675 //plist.FindObject("MTimeEffectiveOnTime")->Clear();
676 }
677 else
678 {
679 // This is the simplest way to remove the two object from the parlist
680 delete halphaoff;
681 delete hfsoff;
682 }
683
684 // ------------- Loop On Data --------------------
685 MReadReports readondata;
686 readondata.AddTree("Events", "MTime.", MReadReports::kMaster);
687 readondata.AddTree("Drive", MReadReports::kRequired);
688 readondata.AddTree("Starguider", MReadReports::kRequired);
689 readondata.AddTree("EffectiveOnTime");
690
691 MReadMarsFile readonmc("Events");
692 readonmc.DisableAutoScheme();
693
694 MRead &readon = fIsMonteCarlo ? (MRead&)readonmc : (MRead&)readondata;
695 set.AddFilesOn(readon);
696
697 scalc.SetMode(MSrcPosCalc::kDefault);
698
699 MFillH fill1b("MHHillasOnPre [MHHillas]", "MHillas", "FillHillasPre");
700 MFillH fill2b("MHHillasOnPost [MHHillas]", "MHillas", "FillHillasPost");
701 MFillH fill3b("MHVsSizeOnPost [MHVsSize]", "MHillasSrc", "FillVsSizePost");
702 MFillH fill4b("MHHilExtOnPost [MHHillasExt]", "MHillasSrc", "FillHilExtPost");
703 MFillH fill5b("MHHilSrcOnPost [MHHillasSrc]", "MHillasSrc", "FillHilSrcPost");
704 MFillH fill6b("MHImgParOnPost [MHImagePar]", "MImagePar", "FillImgParPost");
705 MFillH fill7b("MHNewParOnPost [MHNewImagePar]", "MNewImagePar", "FillNewParPost");
706 fill1b.SetNameTab("PreCut");
707 fill2b.SetNameTab("PostCut");
708 fill3b.SetNameTab("VsSize");
709 fill4b.SetNameTab("HilExt");
710 fill5b.SetNameTab("HilSrc");
711 fill6b.SetNameTab("ImgPar");
712 fill7b.SetNameTab("NewPar");
713 fill1b.SetDrawOption(set.HasOffSequences()||fIsWobble?"same":"");
714 fill2b.SetDrawOption(set.HasOffSequences()||fIsWobble?"same":"");
715 fill3b.SetDrawOption(set.HasOffSequences()||fIsWobble?"same":"");
716 fill4b.SetDrawOption(set.HasOffSequences()||fIsWobble?"same":"");
717 fill5b.SetDrawOption(set.HasOffSequences()||fIsWobble?"same":"");
718 fill6b.SetDrawOption(set.HasOffSequences()||fIsWobble?"same":"");
719 fill7b.SetDrawOption(set.HasOffSequences()||fIsWobble?"same":"");
720
721 /*
722 MHVsTime hvs("MEffectiveOnTime.fVal");
723 hvs.SetTitle("Effective On-Time vs. Time;;T_{on}");
724 MFillH fillvs(&hvs, "MTimeEffectiveOnTime", "FillOnTime");
725 fillvs.SetNameTab("OnTime");
726 */
727 MH3 hvs("MPointingPos.fZd");
728 hvs.SetName("Theta");
729 hvs.SetTitle("Effective On-Time vs. Zenith Angle;\\Theta [\\circ];T_{on} [s]");
730 MFillH fillvs(&hvs, "", "FillOnTime");
731 if (!fIsMonteCarlo)
732 fillvs.SetWeight("MEffectiveOnTime");
733 fillvs.SetNameTab("OnTime");
734
735 /*
736 MParameterD weight;
737 weight.SetVal(-1);
738 fill2a.SetWeight(&weight);
739 fill3a.SetWeight(&weight);
740 fill4a.SetWeight(&weight);
741 fill5a.SetWeight(&weight);
742 fill6a.SetWeight(&weight);
743 fill7a.SetWeight(&weight);
744 if (fSubstraction)
745 {
746 fill2a.SetNameTab("PostCut-");
747 fill3a.SetNameTab("VsSize-");
748 fill4a.SetNameTab("HilExt-");
749 fill5a.SetNameTab("HilSrc-");
750 fill6a.SetNameTab("ImgPar-");
751 fill7a.SetNameTab("NewPar-");
752 }
753 */
754 MHAlpha *halphaon=CreateNewHist(plist);
755 MFillH falpha2(halphaon, "", "FillHist");
756 MH *hfs=CreateNewHistFS(plist);
757 MFillH ffs2(hfs, "MHillas", "FillFS");
758 MFillH fillphi("MHPhi", "", "FillPhi");
759 fillphi.SetDrawOption("anticut");
760
761 tlist.Replace(&readon);
762 if (!fWriteOnly)
763 {
764 tlist2.Replace(&fill1b);
765/* if (fIsWobble)
766 {
767 tlist2.AddToListAfter(&fill2b, &fill2a);
768 tlist2.AddToListAfter(&fill3b, &fill3a);
769 }
770 else
771 */
772 tlist2.Replace(&fill2b);
773 if (fFullDisplay)
774 {
775 tlist2.Replace(&fill3b);
776 tlist2.Replace(&fill4b);
777 tlist2.Replace(&fill5b);
778 tlist2.Replace(&fill6b);
779 tlist2.Replace(&fill7b);
780 }
781 tlist2.Replace(&falpha2);
782 if (!fIsWobble || !fNameHist.IsNull())
783 tlist2.Replace(&ffs2);
784 if (fIsWobble)
785 {
786 tlist2.AddToListAfter(&fillphi, &falpha2);
787 if (!fNameHist.IsNull())
788 tlist2.RemoveFromList(&ffs);
789 }
790
791 if (!fIsMonteCarlo)
792 tlist.AddToList(&fillvs, "EffectiveOnTime");
793 else
794 tlist2.AddToListBefore(&fillvs, &scalc);
795 }
796
797 par.SetVal(1);
798
799 // Execute first analysis
800 if (!evtloop.Eventloop(fMaxEvents))
801 {
802 *fLog << err << GetDescriptor() << ": Processing of on-sequences failed." << endl;
803 return kFALSE;
804 }
805
806 if (write0)
807 delete write0;
808 if (write1)
809 delete write1;
810
811 // FIXME: Perform fit and plot energy dependant alpha plots
812 // and fit result to new tabs!
813 if (!WriteResult(plist, set.GetNumAnalysis()))
814 return kFALSE;
815
816 *fLog << all << GetDescriptor() << ": Done." << endl;
817 *fLog << endl << endl;
818
819 return kTRUE;
820}
Note: See TracBrowser for help on using the repository browser.