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

Last change on this file was 7432, checked in by tbretz, 19 years ago
*** empty log message ***
File size: 25.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// 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// "MHadronness" 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}
368
369// --------------------------------------------------------------------------
370//
371// Create a new instance of an object with name name of class
372// type fNameHist in parlist. It must derive from MHAlpha.
373// Call ForceUsingSize for it and return its pointer.
374// If something fails NULL is returned.
375//
376MHAlpha *MJCut::CreateNewHist(MParList &plist, const char *name) const
377{
378 TClass *cls = gROOT->GetClass(fNameHist);
379 if (!cls)
380 {
381 *fLog << err << "Class " << fNameHist << " not found in dictionary... abort." << endl;
382 return NULL;
383 }
384 if (!cls->InheritsFrom(MHAlpha::Class()))
385 {
386 *fLog << err << "Class " << fNameHist << " doesn't inherit from MHAlpha... abort." << endl;
387 return NULL;
388 }
389
390 const TString objname(Form("Hist%s", name));
391 MHAlpha *h = (MHAlpha*)plist.FindCreateObj(fNameHist, objname);
392 if (!h)
393 return NULL;
394
395 h->ForceUsingSize();
396
397 return h;
398}
399
400// --------------------------------------------------------------------------
401//
402// Create a new instance of an object with name name of class
403// type fNameHistFS in parlist. It must derive from MHFalseSource
404// If something fails NULL is returned.
405//
406MH *MJCut::CreateNewHistFS(MParList &plist, const char *name) const
407{
408 const TString cname(fNameHistFS.IsNull()?"MHFalseSource":fNameHistFS);
409
410 TClass *cls = gROOT->GetClass(cname);
411 if (!cls)
412 {
413 *fLog << err << "Class " << cname << " not found in dictionary... abort." << endl;
414 return NULL;
415 }
416 if (!cls->InheritsFrom("MHFalseSource"))
417 {
418 *fLog << err << "Class " << cname << " doesn't inherit from MHFalseSource... abort." << endl;
419 return NULL;
420 }
421
422 const TString objname(Form("FS%s", name));
423 return (MH*)plist.FindCreateObj(cname, objname);
424}
425
426Bool_t MJCut::Process(const MDataSet &set)
427{
428 if (!set.IsValid())
429 {
430 *fLog << err << "ERROR - DataSet invalid!" << endl;
431 return kFALSE;
432 }
433
434 CheckEnv();
435
436 // --------------------------------------------------------------------------------
437
438 *fLog << inf;
439 fLog->Separator(GetDescriptor());
440 *fLog << "Perform cuts for data set " << set.GetName() << endl;
441 *fLog << endl;
442
443 // --------------------------------------------------------------------------------
444
445 // Setup Parlist
446 MParList plist;
447 plist.AddToList(this); // take care of fDisplay!
448
449 MParameterI par("DataType");
450 plist.AddToList(&par);
451
452 // Setup Tasklist
453 MTaskList tlist;
454 plist.AddToList(&tlist);
455
456 // La Palma Magic1
457 MObservatory obs;
458 plist.AddToList(&obs);
459
460 // Possible source position (eg. Wobble Mode)
461 MPointingPos source("MSourcePos");
462 if (set.HasSource())
463 {
464 if (!set.GetSourcePos(source))
465 return kFALSE;
466 plist.AddToList(&source);
467 *fLog << all << "Using Source Position: ";
468 source.Print("RaDec");
469 }
470 else
471 *fLog << inf << "No source position applied..." << endl;
472
473 // Initialize default binnings
474 MBinning bins1(18, 0, 90, "BinningAlpha", "lin");
475 MBinning bins2(15, 10, 1e6 , "BinningSize", "log");
476 MBinning bins3(67, -0.005, 0.665, "BinningTheta", "asin");
477 MBinning bins4("BinningFalseSource");
478 MBinning bins5("BinningWidth");
479 MBinning bins6("BinningLength");
480 MBinning bins7("BinningDist");
481 MBinning bins8("BinningMaxDist");
482 MBinning bins9("BinningM3Long");
483 MBinning bins0("BinningConc1");
484 plist.AddToList(&bins1);
485 plist.AddToList(&bins2);
486 plist.AddToList(&bins3);
487 plist.AddToList(&bins4);
488 plist.AddToList(&bins5);
489 plist.AddToList(&bins6);
490 plist.AddToList(&bins7);
491 plist.AddToList(&bins8);
492 plist.AddToList(&bins9);
493 plist.AddToList(&bins0);
494 //plist.AddToList(&binsa);
495
496 // --------------------------------------------------------------------------------
497
498 // Setup fitter and histograms
499 MAlphaFitter fit;
500 plist.AddToList(&fit);
501 if (fIsWobble)
502 fit.SetScaleMode(MAlphaFitter::kNone);
503
504 MHAlpha *halphaoff = CreateNewHist(plist, "Off");
505 MFillH falpha(halphaoff, "", "FillHist");
506 MH *hfsoff = CreateNewHistFS(plist, "Off");
507 MFillH ffs(hfsoff, "MHillas", "FillFS");
508
509 // FIXME: If fPathIn read cuts and energy estimator from file!
510 MContinue cont0("", "Cut0");
511 MContinue cont1("", "Cut1");
512 MContinue cont2("", "Cut2");
513 MContinue cont3("", "Cut3");
514 cont0.SetAllowEmpty();
515 cont1.SetAllowEmpty();
516 cont2.SetAllowEmpty();
517 cont3.SetAllowEmpty();
518
519 // ------------- Loop Off Data --------------------
520 MReadReports readoffdata;
521 readoffdata.AddTree("Events", "MTime.", MReadReports::kMaster);
522 readoffdata.AddTree("Drive", MReadReports::kRequired);
523 readoffdata.AddTree("Starguider", MReadReports::kRequired);
524 readoffdata.AddTree("EffectiveOnTime");
525
526 MReadMarsFile readoffmc("Events");
527 readoffmc.DisableAutoScheme();
528
529 MRead &readoff = fIsMonteCarlo ? (MRead&)readoffmc : (MRead&)readoffdata;
530 if (fIsWobble)
531 set.AddFilesOn(readoff);
532 else
533 set.AddFilesOff(readoff);
534
535 const TString path(Form("%s/", fPathOut.Data()));
536 TString fname0(path);
537 TString fname1(path);
538 fname0 += fNameSummary.IsNull() ? (TString) Form("ganymed%08d-summary.root", set.GetNumAnalysis()) : fNameSummary;
539 fname1 += fNameResult.IsNull() ? (TString) Form("ganymed%08d.root", set.GetNumAnalysis()) : fNameResult;
540
541 MWriteRootFile *write0 = CanStoreSummary() ? new MWriteRootFile(fPathOut.IsNull()?0:fname0.Data(), fOverwrite?"RECREATE":"NEW") : 0;
542 MWriteRootFile *write1 = CanStoreResult() ? new MWriteRootFile(fPathOut.IsNull()?0:fname1.Data(), fOverwrite?"RECREATE":"NEW") : 0;
543 SetupWriter(write0, "WriteAfterCut0");
544 SetupWriter(write1, "WriteAfterCut3");
545
546/*
547 MEnergyEstimate est;
548
549 MTaskEnv taskenv1("EstimateEnergy");
550 taskenv1.SetDefault(fEstimateEnergy ? fEstimateEnergy : &est);
551 */
552 MTaskEnv taskenv2("CalcHadronness");
553 taskenv2.SetDefault(fCalcHadronness);
554
555 MTaskEnv taskenv3("CalcDisp");
556 taskenv3.SetDefault(fCalcDisp);
557
558 MFillH fill1a("MHHillasOffPre [MHHillas]", "MHillas", "FillHillasPre");
559 MFillH fill2a("MHHillasOffPost [MHHillas]", "MHillas", "FillHillasPost");
560 MFillH fill3a("MHVsSizeOffPost [MHVsSize]", "MHillasSrc", "FillVsSizePost");
561 MFillH fill4a("MHHilExtOffPost [MHHillasExt]", "MHillasSrc", "FillHilExtPost");
562 MFillH fill5a("MHHilSrcOffPost [MHHillasSrc]", "MHillasSrc", "FillHilSrcPost");
563 MFillH fill6a("MHImgParOffPost [MHImagePar]", "MImagePar", "FillImgParPost");
564 MFillH fill7a("MHNewParOffPost [MHNewImagePar]", "MNewImagePar", "FillNewParPost");
565 fill1a.SetNameTab("PreCut");
566 fill2a.SetNameTab("PostCut");
567 fill3a.SetNameTab("VsSize");
568 fill4a.SetNameTab("HilExt");
569 fill5a.SetNameTab("HilSrc");
570 fill6a.SetNameTab("ImgPar");
571 fill7a.SetNameTab("NewPar");
572
573 MPrint print2("MEffectiveOnTime");
574 print2.EnableSkip();
575
576 // How to get source position from off- and on-data?
577 MSrcPosCalc scalc;
578 scalc.SetMode(fIsWobble?MSrcPosCalc::kWobble:MSrcPosCalc::kOffData); /********************/
579
580 MSrcPosCorrect scor;
581
582 MHillasCalc hcalc;
583 MHillasCalc hcalc2("MHillasCalcAnti");
584 hcalc.SetFlags(MHillasCalc::kCalcHillasSrc);
585 hcalc2.SetFlags(MHillasCalc::kCalcHillasSrc);
586 hcalc2.SetNameHillasSrc("MHillasSrcAnti");
587 hcalc2.SetNameSrcPosCam("MSrcPosAnti");
588
589 MTaskList tlist2;
590 tlist2.AddToList(&scalc);
591 tlist2.AddToList(&scor);
592 tlist2.AddToList(&hcalc);
593 if (fIsWobble)
594 tlist2.AddToList(&hcalc2);
595 //tlist2.AddToList(&taskenv1);
596 tlist2.AddToList(&taskenv2);
597 tlist2.AddToList(&taskenv3);
598 tlist2.AddToList(&cont0);
599 if (write0)
600 tlist2.AddToList(write0);
601 if (!fWriteOnly)
602 tlist2.AddToList(&fill1a);
603 tlist2.AddToList(&cont1);
604 if (!fWriteOnly && (!fIsWobble || !fNameHistFS.IsNull()))
605 tlist2.AddToList(&ffs);
606 tlist2.AddToList(&cont2);
607 if (!fWriteOnly)
608 {
609 tlist2.AddToList(&fill2a);
610 if (fFullDisplay)
611 {
612 tlist2.AddToList(&fill3a);
613 tlist2.AddToList(&fill4a);
614 tlist2.AddToList(&fill5a);
615 tlist2.AddToList(&fill6a);
616 tlist2.AddToList(&fill7a);
617 }
618 }
619 if (!fWriteOnly)
620 tlist2.AddToList(&falpha);
621 tlist2.AddToList(&cont3);
622 if (write1)
623 tlist2.AddToList(write1);
624
625 MPointingDevCalc devcalc;
626
627 tlist.AddToList(&readoff);
628 if (gLog.GetDebugLevel()>4)
629 tlist.AddToList(&print2, "EffectiveOnTime");
630 tlist.AddToList(&devcalc, "Starguider");
631 tlist.AddToList(&tlist2, "Events");
632
633 par.SetVal(0);
634
635 // Create and setup the eventloop
636 MEvtLoop evtloop(fName);
637 evtloop.SetParList(&plist);
638 evtloop.SetDisplay(fDisplay);
639 evtloop.SetLogStream(fLog);
640 if (!SetupEnv(evtloop))
641 return kFALSE;
642
643 TObjArray cont;
644 cont.Add(&cont0);
645 cont.Add(&cont1);
646 cont.Add(&cont2);
647 cont.Add(&cont3);
648 //if (taskenv1.GetTask())
649 // cont.Add(taskenv1.GetTask());
650 if (taskenv2.GetTask())
651 cont.Add(taskenv2.GetTask());
652 if (taskenv3.GetTask())
653 cont.Add(taskenv3.GetTask());
654
655 if (!WriteTasks(set.GetNumAnalysis(), cont))
656 return kFALSE;
657
658 if (set.HasOffSequences() || fIsWobble)
659 {
660 // Execute first analysis
661 if (!evtloop.Eventloop(fMaxEvents))
662 {
663 *fLog << err << GetDescriptor() << ": Processing of off-sequences failed." << endl;
664 return kFALSE;
665 }
666
667 if (!evtloop.GetDisplay())
668 {
669 *fLog << err << GetDescriptor() << ": Execution stopped by user." << endl;
670 return kFALSE;
671 }
672
673 //plist.FindObject("MTimeEffectiveOnTime")->Clear();
674 }
675 else
676 {
677 // This is the simplest way to remove the two object from the parlist
678 delete halphaoff;
679 delete hfsoff;
680 }
681
682 // ------------- Loop On Data --------------------
683 MReadReports readondata;
684 readondata.AddTree("Events", "MTime.", MReadReports::kMaster);
685 readondata.AddTree("Drive", MReadReports::kRequired);
686 readondata.AddTree("Starguider", MReadReports::kRequired);
687 readondata.AddTree("EffectiveOnTime");
688
689 MReadMarsFile readonmc("Events");
690 readonmc.DisableAutoScheme();
691
692 MRead &readon = fIsMonteCarlo ? (MRead&)readonmc : (MRead&)readondata;
693 set.AddFilesOn(readon);
694
695 scalc.SetMode(MSrcPosCalc::kDefault);
696
697 MFillH fill1b("MHHillasOnPre [MHHillas]", "MHillas", "FillHillasPre");
698 MFillH fill2b("MHHillasOnPost [MHHillas]", "MHillas", "FillHillasPost");
699 MFillH fill3b("MHVsSizeOnPost [MHVsSize]", "MHillasSrc", "FillVsSizePost");
700 MFillH fill4b("MHHilExtOnPost [MHHillasExt]", "MHillasSrc", "FillHilExtPost");
701 MFillH fill5b("MHHilSrcOnPost [MHHillasSrc]", "MHillasSrc", "FillHilSrcPost");
702 MFillH fill6b("MHImgParOnPost [MHImagePar]", "MImagePar", "FillImgParPost");
703 MFillH fill7b("MHNewParOnPost [MHNewImagePar]", "MNewImagePar", "FillNewParPost");
704 fill1b.SetNameTab("PreCut");
705 fill2b.SetNameTab("PostCut");
706 fill3b.SetNameTab("VsSize");
707 fill4b.SetNameTab("HilExt");
708 fill5b.SetNameTab("HilSrc");
709 fill6b.SetNameTab("ImgPar");
710 fill7b.SetNameTab("NewPar");
711 fill1b.SetDrawOption(set.HasOffSequences()||fIsWobble?"same":"");
712 fill2b.SetDrawOption(set.HasOffSequences()||fIsWobble?"same":"");
713 fill3b.SetDrawOption(set.HasOffSequences()||fIsWobble?"same":"");
714 fill4b.SetDrawOption(set.HasOffSequences()||fIsWobble?"same":"");
715 fill5b.SetDrawOption(set.HasOffSequences()||fIsWobble?"same":"");
716 fill6b.SetDrawOption(set.HasOffSequences()||fIsWobble?"same":"");
717 fill7b.SetDrawOption(set.HasOffSequences()||fIsWobble?"same":"");
718
719 /*
720 MHVsTime hvs("MEffectiveOnTime.fVal");
721 hvs.SetTitle("Effective On-Time vs. Time;;T_{on}");
722 MFillH fillvs(&hvs, "MTimeEffectiveOnTime", "FillOnTime");
723 fillvs.SetNameTab("OnTime");
724 */
725 MH3 hvs("MPointingPos.fZd");
726 hvs.SetName("Theta");
727 hvs.SetTitle("Effective On-Time vs. Zenith Angle;\\Theta [\\circ];T_{on} [s]");
728 MFillH fillvs(&hvs, "", "FillOnTime");
729 if (!fIsMonteCarlo)
730 fillvs.SetWeight("MEffectiveOnTime");
731 fillvs.SetNameTab("OnTime");
732
733 /*
734 MParameterD weight;
735 weight.SetVal(-1);
736 fill2a.SetWeight(&weight);
737 fill3a.SetWeight(&weight);
738 fill4a.SetWeight(&weight);
739 fill5a.SetWeight(&weight);
740 fill6a.SetWeight(&weight);
741 fill7a.SetWeight(&weight);
742 if (fSubstraction)
743 {
744 fill2a.SetNameTab("PostCut-");
745 fill3a.SetNameTab("VsSize-");
746 fill4a.SetNameTab("HilExt-");
747 fill5a.SetNameTab("HilSrc-");
748 fill6a.SetNameTab("ImgPar-");
749 fill7a.SetNameTab("NewPar-");
750 }
751 */
752 MHAlpha *halphaon=CreateNewHist(plist);
753 MFillH falpha2(halphaon, "", "FillHist");
754 MH *hfs=CreateNewHistFS(plist);
755 MFillH ffs2(hfs, "MHillas", "FillFS");
756 MFillH fillphi("MHPhi", "", "FillPhi");
757 fillphi.SetDrawOption("anticut");
758
759 tlist.Replace(&readon);
760 if (!fWriteOnly)
761 {
762 tlist2.Replace(&fill1b);
763/* if (fIsWobble)
764 {
765 tlist2.AddToListAfter(&fill2b, &fill2a);
766 tlist2.AddToListAfter(&fill3b, &fill3a);
767 }
768 else
769 */
770 tlist2.Replace(&fill2b);
771 if (fFullDisplay)
772 {
773 tlist2.Replace(&fill3b);
774 tlist2.Replace(&fill4b);
775 tlist2.Replace(&fill5b);
776 tlist2.Replace(&fill6b);
777 tlist2.Replace(&fill7b);
778 }
779 tlist2.Replace(&falpha2);
780 if (!fIsWobble || !fNameHist.IsNull())
781 tlist2.Replace(&ffs2);
782 if (fIsWobble)
783 {
784 tlist2.AddToListAfter(&fillphi, &falpha2);
785 if (!fNameHist.IsNull())
786 tlist2.RemoveFromList(&ffs);
787 }
788
789 if (!fIsMonteCarlo)
790 tlist.AddToList(&fillvs, "EffectiveOnTime");
791 else
792 tlist2.AddToListBefore(&fillvs, &scalc);
793 }
794
795 par.SetVal(1);
796
797 // Execute first analysis
798 if (!evtloop.Eventloop(fMaxEvents))
799 {
800 *fLog << err << GetDescriptor() << ": Processing of on-sequences failed." << endl;
801 return kFALSE;
802 }
803
804 if (write0)
805 delete write0;
806 if (write1)
807 delete write1;
808
809 // FIXME: Perform fit and plot energy dependant alpha plots
810 // and fit result to new tabs!
811 if (!WriteResult(plist, set.GetNumAnalysis()))
812 return kFALSE;
813
814 *fLog << all << GetDescriptor() << ": Done." << endl;
815 *fLog << endl << endl;
816
817 return kTRUE;
818}
Note: See TracBrowser for help on using the repository browser.