source: branches/Mars_McMismatchStudy/mjobs/MJCut.cc@ 18788

Last change on this file since 18788 was 14893, checked in by tbretz, 12 years ago
Adapted the call to MHSrcPosCam
File size: 34.1 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-2008
21!
22!
23\* ======================================================================== */
24
25/////////////////////////////////////////////////////////////////////////////
26//
27// MJCut
28//
29/////////////////////////////////////////////////////////////////////////////
30#include "MJCut.h"
31
32// Root
33#include <TEnv.h>
34#include <TFile.h>
35#include <TClass.h>
36
37// Environment
38#include "MLog.h"
39#include "MLogManip.h"
40
41// Eventloop
42#include "MParList.h"
43#include "MTaskList.h"
44#include "MEvtLoop.h"
45
46// Display
47#include "MStatusDisplay.h"
48
49// Tasks
50#include "MReadReports.h"
51#include "MReadMarsFile.h"
52#include "MPrint.h"
53#include "MContinue.h"
54#include "MTaskEnv.h"
55#include "MPointingDevCalc.h"
56#include "MParameterCalc.h"
57#include "MSrcPosRndm.h"
58#include "MSrcPosCalc.h"
59#include "MSrcPosCorrect.h"
60#include "MHillasCalc.h"
61#include "MFillH.h"
62#include "MWriteRootFile.h"
63
64// Filter
65#include "MFDataPhrase.h"
66
67// Fit signal environment
68#include "../mhflux/MAlphaFitter.h"
69#include "../mhflux/MHAlpha.h"
70
71// Containers
72#include "MH3.h"
73#include "MBinning.h"
74#include "MDataSet.h"
75#include "MParameters.h"
76#include "MPointingPos.h"
77#include "MObservatory.h"
78#include "MHSrcPosCam.h"
79
80ClassImp(MJCut);
81
82using namespace std;
83
84// --------------------------------------------------------------------------
85//
86// Default constructor. Set defaults for fStoreSummary, fStoreresult,
87// fWriteOnly, fFullDisplay to kFALSE and initialize
88// /*fEstimateEnergy and*/ fCalcHadronness with NULL.
89//
90MJCut::MJCut(const char *name, const char *title)
91 : fStoreSummary(kFALSE), fStoreResult(kTRUE), fWriteOnly(kFALSE),
92 fFullDisplay(kTRUE), fRndmSrcPos(kFALSE), fNumOffSourcePos(3),
93 fNameHist("MHThetaSq"),
94 fCalcHadronness(0), fCalcDisp(0), fEstimateEnergy(0)
95{
96 fName = name ? name : "MJCut";
97 fTitle = title ? title : "Standard program to perform g/h-separation cuts";
98}
99
100// --------------------------------------------------------------------------
101//
102// Destructor. Delete fEstimateEnergy and fCalcHadronness if != NULL
103//
104MJCut::~MJCut()
105{
106 if (fEstimateEnergy)
107 delete fEstimateEnergy;
108 if (fCalcHadronness)
109 delete fCalcHadronness;
110 if (fCalcDisp)
111 delete fCalcDisp;
112}
113
114// --------------------------------------------------------------------------
115//
116// Set the name of the summary file (events after cut0)
117// If you give a name the storage of this file is enabled implicitly.
118// If you give no filename the storage is neither enabled nor disabled,
119// but the storage file name is reset.
120// If no filename is set the default filename is used.
121// You can explicitly enable or disable the storage using EnableStoreOf*()
122// The default argument is no filename.
123//
124void MJCut::SetNameSummaryFile(const char *name)
125{
126 fNameSummary=name;
127 if (!fNameSummary.IsNull())
128 EnableStorageOfSummary();
129}
130
131// --------------------------------------------------------------------------
132//
133// Set the name of the summary file (events after cut3)
134// If you give a name the storage of this file is enabled implicitly.
135// If you give no filename the storage is neither enabled nor disabled,
136// but the storage file name is reset.
137// If no filename is set the default filename is used.
138// You can explicitly enable or disable the storage using EnableStoreOf*()
139// The default argument is no filename.
140//
141void MJCut::SetNameResultFile(const char *name)
142{
143 fNameResult=name;
144 if (!fNameResult.IsNull())
145 EnableStorageOfResult();
146}
147
148// --------------------------------------------------------------------------
149//
150// Setup a task estimating the energy. The given task is cloned.
151//
152/*
153void MJCut::SetEnergyEstimator(const MTask *task)
154{
155 if (fEstimateEnergy)
156 delete fEstimateEnergy;
157 fEstimateEnergy = task ? (MTask*)task->Clone() : 0;
158}
159*/
160
161// --------------------------------------------------------------------------
162//
163// Setup a task calculating the hadronness. The given task is cloned.
164//
165void MJCut::SetHadronnessCalculator(const MTask *task)
166{
167 if (fCalcHadronness)
168 delete fCalcHadronness;
169 fCalcHadronness = task ? (MTask*)task->Clone() : 0;
170}
171
172// --------------------------------------------------------------------------
173//
174// Setup a task calculating disp. The given task is cloned.
175//
176void MJCut::SetDispCalculator(const MTask *task)
177{
178 if (fCalcDisp)
179 delete fCalcDisp;
180 fCalcDisp = task ? (MTask*)task->Clone() : 0;
181}
182
183// --------------------------------------------------------------------------
184//
185// Setup a task estimating the eneryg. The given task is cloned.
186//
187void MJCut::SetEnergyEstimator(const MTask *task)
188{
189 if (fEstimateEnergy)
190 delete fEstimateEnergy;
191 fEstimateEnergy = task ? (MTask*)task->Clone() : 0;
192}
193
194// --------------------------------------------------------------------------
195//
196// return fOutputPath+"/ganymed%08d.root", num
197//
198TString MJCut::GetOutputFile(UInt_t num) const
199{
200 TString p(fPathOut);
201 p += "/";
202 p += fNameOutput.IsNull() ? Form("ganymed%08d.root", num) : fNameOutput.Data();
203 gSystem->ExpandPathName(p);
204 return p;
205}
206
207/*
208Bool_t MJCut::ReadTasks(const char *fname, MTask* &env1, MTask* &env2) const
209{
210 // const TString fname = Form("%s/calib%08d.root", fPathIn.Data(), fSequence.GetSequence());
211
212 *fLog << inf << "Reading from file: " << fname << endl;
213
214 TFile file(fname, "READ");
215 if (!file.IsOpen())
216 {
217 *fLog << err << dbginf << "ERROR - Could not open file " << fname << endl;
218 return kFALSE;
219 }
220
221 TObject *o = file.Get("EstimateEnergy");
222 if (o && !o->InheritsFrom(MTask::Class()))
223 {
224 *fLog << err << dbginf << "ERROR - EstimateEnergy read from " << fname << " doesn't inherit from MTask!" << endl;
225 return kFALSE;
226 }
227 env1 = o ? (MTask*)o->Clone() : FNULL;
228
229 o = file.Get("CalcHadronness");
230 if (o && !o->InheritsFrom(MTask::Class()))
231 {
232 *fLog << err << dbginf << "ERROR - CalcHadronness read from " << fname << " doesn't inherit from MTask!" << endl;
233 return kFALSE;
234 }
235 env2 = o ? (MTask*)o->Clone() : NULL;
236
237 return kTRUE;
238}
239*/
240
241// --------------------------------------------------------------------------
242//
243// Write the tasks in cont to the file corresponding to analysis number num,
244// see GetOutputFile()
245//
246Bool_t MJCut::WriteTasks(UInt_t num, TObjArray &cont) const
247{
248 if (fPathOut.IsNull())
249 {
250 *fLog << inf << "No output path specified via SetPathOut - no output written." << endl;
251 return kTRUE;
252 }
253
254 const TString oname(GetOutputFile(num));
255
256 *fLog << inf << "Writing to file: " << oname << endl;
257
258 TFile *file = 0;
259 if (fNameResult.IsNull() && fStoreResult)
260 {
261 // If no special name fNameResult is given but the result should be stored,
262 // search for a file named oname (fNameOutput or ganymed%08d.root)
263 file = (TFile*)gROOT->GetListOfFiles()->FindObject(oname);
264 if (!file)
265 {
266 *fLog << err << "ERROR - File " << oname << " expected to be open, but isn't..." << endl;
267 return kFALSE;
268 }
269
270 file->cd();
271 }
272 else
273 {
274 // If special name fNameResult is given or no storage requested,
275 // open a new file named oname (fNameOutput or ganymed%08d.root)
276 file = TFile::Open(oname, fOverwrite?"RECREATE":"NEW", "File created by MJCut", 9);
277
278 if (!file)
279 {
280 *fLog << err << "ERROR - Couldn't open new file " << oname << " for writing..." << endl;
281 return kFALSE;
282 }
283 }
284
285 // Write container to output file
286 const Bool_t rc = WriteContainer(cont);
287
288 // If file was newly opened, close it
289 if (!(fNameResult.IsNull() && fStoreResult))
290 delete file;
291
292 // Return result of writing
293 return rc;
294}
295
296// --------------------------------------------------------------------------
297//
298// Write the result plots and other results to the file corresponding to
299// analysis number num, see GetOutputFile()
300//
301Bool_t MJCut::WriteResult(const MParList &plist, const MDataSet &set) const
302{
303 TObjArray arr;
304
305 // Save all MBinnings
306 TIter Next(plist);
307 TObject *o=0;
308 while ((o=Next()))
309 if (o->InheritsFrom(MBinning::Class()))
310 arr.Add(o);
311
312 // Save also the result, not only the setup
313 const MHAlpha *halpha = (MHAlpha*)plist.FindObject("Hist", "MHAlpha");
314 if (halpha)
315 arr.Add((TObject*)(&halpha->GetAlphaFitter()));
316
317 // Save also the dataset, environment and display
318 arr.Add(const_cast<MDataSet*>(&set));
319 arr.Add(const_cast<TEnv*>(GetEnv()));
320 if (fDisplay)
321 arr.Add(fDisplay);
322
323 TNamed cmdline("CommandLine", fCommandLine.Data());
324 arr.Add(&cmdline);
325
326 const Int_t num = set.GetNumAnalysis();
327 const TString fname(fNameOutput.IsNull() ? Form("ganymed%08d.root", num) : fNameOutput.Data());
328
329 // If requested, write to already open output file
330 if (fNameResult.IsNull() && fStoreResult)
331 {
332 *fLog << inf2 << "Searching for open file " << fname << "... " << flush;
333
334 TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject(fname);
335 if (file)
336 {
337 *fLog << "found." << endl;
338 file->cd();
339 return WriteContainer(arr);
340 }
341 *fLog << "not found." << endl;
342 }
343
344 return WriteContainer(arr, fname, "UPDATE");
345}
346
347// --------------------------------------------------------------------------
348//
349// MJCut allows to setup several option by a resource file:
350// MJCut.WriteSummary: yes, no
351// MJCut.SummaryFile: filename
352// MJCut.WriteResult: yes, no
353// MJCut.ResultFile: filename
354// MJCut.HistName: MHAlpha
355//
356Bool_t MJCut::CheckEnvLocal()
357{
358 const TString f0(GetEnv("SummaryFile", ""));
359 const TString f1(GetEnv("ResultFile", ""));
360 if (!f0.IsNull())
361 SetNameSummaryFile(f0);
362 if (!f1.IsNull())
363 SetNameResultFile(f1);
364
365 EnableStorageOfSummary(GetEnv("SummaryFile", fStoreSummary));
366 EnableStorageOfResult(GetEnv("ResultFile", fStoreResult));
367 EnableFullDisplay(GetEnv("FullDisplay", fFullDisplay));
368 EnableRandomSrcPos(GetEnv("RandomSourcePosition", fRndmSrcPos));
369
370 fNumOffSourcePos = GetEnv("NumOffSourcePositions", (Int_t)fNumOffSourcePos);
371 //EnableSubstraction(GetEnv("HistogramSubstraction", fSubstraction));
372
373 SetNameHist(GetEnv("NameHist", fNameHist));
374 SetNameHistFS(GetEnv("NameHistFS", fNameHistFS));
375
376 return kTRUE;
377}
378
379// --------------------------------------------------------------------------
380//
381// Setup write to write:
382// container tree optional?
383// -------------- ---------- -----------
384// "MHillas" to "Events"
385// "MHillasSrc" to "Events"
386// "Hadronness" to "Events" yes
387// "MEnergyEst" to "Events" yes
388// "DataType" to "Events"
389//
390void MJCut::SetupWriter(MWriteRootFile *write, const char *name) const
391{
392 if (!write)
393 return;
394
395 write->SetName(name);
396 write->AddContainer("MHillas", "Events");
397 write->AddContainer("MHillasSrc", "Events");
398 write->AddContainer("MHillasExt", "Events");
399 write->AddContainer("MPointingPos", "Events");
400 write->AddContainer("MHillasSrcAnti", "Events", kFALSE);
401 write->AddContainer("MImagePar", "Events", kFALSE);
402 write->AddContainer("MNewImagePar", "Events", kFALSE);
403 write->AddContainer("MNewImagePar2", "Events", kFALSE);
404 write->AddContainer("Hadronness", "Events", kFALSE);
405 write->AddContainer("MSrcPosCam", "Events", kFALSE);
406 write->AddContainer("MSrcPosAnti", "Events", kFALSE);
407 write->AddContainer("ThetaSquared", "Events", kFALSE);
408 write->AddContainer("OpticalAxis", "Events", kFALSE);
409 write->AddContainer("Disp", "Events", kFALSE);
410 write->AddContainer("Ghostbuster", "Events", kFALSE);
411 write->AddContainer("MEnergyEst", "Events", kFALSE);
412 write->AddContainer("MTime", "Events", kFALSE);
413 write->AddContainer("MMcEvt", "Events", kFALSE);
414 write->AddContainer("DataType", "Events");
415 write->AddContainer("FileId", "Events");
416 write->AddContainer("EvtNumber", "Events");
417 // write->AddContainer("MMuonSearchPar", "Events", kFALSE);
418 // write->AddContainer("MMuonCalibPar", "Events", kFALSE);
419}
420
421// --------------------------------------------------------------------------
422//
423// Create a new instance of an object with name name of class
424// type fNameHist in parlist. It must derive from MHAlpha.
425// Call ForceUsingSize for it and return its pointer.
426// If something fails NULL is returned.
427//
428MHAlpha *MJCut::CreateNewHist(MParList &plist, const char *name) const
429{
430 TClass *cls = gROOT->GetClass(fNameHist);
431 if (!cls)
432 {
433 *fLog << err << "Class " << fNameHist << " not found in dictionary... abort." << endl;
434 return NULL;
435 }
436 if (!cls->InheritsFrom(MHAlpha::Class()))
437 {
438 *fLog << err << "Class " << fNameHist << " doesn't inherit from MHAlpha... abort." << endl;
439 return NULL;
440 }
441
442 const TString objname(Form("Hist%s", name));
443 MHAlpha *h = (MHAlpha*)plist.FindCreateObj(fNameHist, objname);
444 if (!h)
445 return NULL;
446
447 h->ForceUsingSize();
448
449 return h;
450}
451
452// --------------------------------------------------------------------------
453//
454// Create a new instance of an object with name name of class
455// type fNameHistFS in parlist. It must derive from MHFalseSource
456// If something fails NULL is returned.
457//
458MH *MJCut::CreateNewHistFS(MParList &plist, const char *name) const
459{
460 const TString cname(fNameHistFS.IsNull()?"MHFalseSource":fNameHistFS.Data());
461
462 TClass *cls = gROOT->GetClass(cname);
463 if (!cls)
464 {
465 *fLog << err << "Class " << cname << " not found in dictionary... abort." << endl;
466 return NULL;
467 }
468 if (!cls->InheritsFrom("MHFalseSource"))
469 {
470 *fLog << err << "Class " << cname << " doesn't inherit from MHFalseSource... abort." << endl;
471 return NULL;
472 }
473
474 const TString objname(Form("FS%s", name));
475 return (MH*)plist.FindCreateObj(cname, objname);
476}
477
478Bool_t MJCut::FillSrcPosCam(const MDataSet &set, MPointingPos &source, MHSrcPosCam &hsrcpos)
479{
480 *fLog << inf;
481 fLog->Separator(GetDescriptor());
482 *fLog << "Filling MHSrcPosCam " << set.GetBaseName() << endl;
483 *fLog << endl;
484
485 // --------------------------------------------------------------------------------
486
487 // Setup Parlist
488 MParList plist;
489 plist.AddToList(this); // take care of fDisplay!
490
491 // Setup Tasklist
492 MTaskList tlist;
493 plist.AddToList(&tlist);
494
495 // La Palma Magic1, Possible source position
496 MObservatory obs;
497 plist.AddToList(&obs);
498 plist.AddToList(&source);
499
500 // Initialize default binnings
501 // MBinning bins1(18, 0, 90, "BinningSrcPosCam", "lin");
502 // plist.AddToList(&bins1);
503
504 // ------------- Loop Off Data --------------------
505 MReadReports read;
506
507 read.EnableAutoScheme();
508 read.AddToBranchList("MTimeEffectiveOnTime.*");
509 read.AddToBranchList("MEffectiveOnTime.*");
510
511 read.AddTree("Events", "MTime.", MReadReports::kMaster);
512 read.AddTree("Drive", MReadReports::kRequired);
513 read.AddTree("Starguider", MReadReports::kRequired);
514 read.AddTree("EffectiveOnTime");
515
516 if (!set.AddFilesOn(read))
517 return kFALSE;
518
519 MFillH fill(&hsrcpos, "MSrcPosCam", "FillSrcPosCam");
520 fill.SetNameTab("SrcPos");
521
522 // How to get source position from off- and on-data?
523 MSrcPosCorrect scor;
524 MSrcPosCalc scalc;
525 scalc.SetMode(MSrcPosCalc::kDefault);
526
527 MPointingDevCalc devcalc;
528
529 tlist.AddToList(&read);
530 tlist.AddToList(&devcalc, "Starguider");
531 tlist.AddToList(&scalc, "Events");
532 tlist.AddToList(&scor, "Events");
533 tlist.AddToList(&fill, "Events");
534
535 // by setting it here it is distributed to all consecutive tasks
536 tlist.SetAccelerator(MTask::kAccDontReset|MTask::kAccDontTime);
537
538 // Create and setup the eventloop
539 MEvtLoop evtloop(fName);
540 evtloop.SetParList(&plist);
541 evtloop.SetDisplay(fDisplay);
542 evtloop.SetLogStream(fLog);
543 if (!SetupEnv(evtloop))
544 return kFALSE;
545
546 // Execute first analysis
547 if (!evtloop.Eventloop(fMaxEvents))
548 {
549 *fLog << err << GetDescriptor() << ": Processing of on-sequences to fill SrcPosCam failed." << endl;
550 return kFALSE;
551 }
552
553 if (!evtloop.GetDisplay())
554 {
555 *fLog << err << GetDescriptor() << ": Execution stopped by user." << endl;
556 return kFALSE;
557 }
558
559 *fLog << all << GetDescriptor() << ": Done." << endl;
560 *fLog << endl << endl;
561
562 return kTRUE;
563}
564
565Int_t MJCut::Process(const MDataSet &set)
566{
567 if (!set.IsValid())
568 {
569 *fLog << err << "ERROR - DataSet invalid!" << endl;
570 return kFALSE;
571 }
572
573 if (!HasWritePermission(GetOutputFile(set.GetNumAnalysis())))
574 return kFALSE;
575
576 if (!CheckEnv())
577 return kFALSE;
578
579 // --------------------------------------------------------------------------------
580
581 if (!set.IsWobbleMode() && fNumOffSourcePos!=1)
582 {
583 *fLog << inf << "No wobble mode, but NumOffSoucePos!=1 (" << fNumOffSourcePos << ")... reset to 1." << endl;
584 fNumOffSourcePos = 1;
585 }
586
587 // Possible source position (eg. Wobble Mode)
588 MPointingPos source("MSourcePos");
589 if (set.HasSource())
590 {
591 if (!set.GetSourcePos(source))
592 return -1;
593 *fLog << all;
594 source.Print("RaDec");
595 }
596 else
597 *fLog << all << "No source position applied..." << endl;
598
599 // ADD A CHECK FOR THE SOURCE POSITION by COMPARISON WITH RA/DEC
600
601 MParList plist;
602
603 MHSrcPosCam hsrcpos(set.IsWobbleMode()?0.4:0);
604 if (!set.IsWobbleMode() && source.IsInitialized() && fRndmSrcPos)
605 {
606 if (!FillSrcPosCam(set, source, hsrcpos))
607 return -2;
608 plist.AddToList(&hsrcpos);
609 }
610
611 // --------------------------------------------------------------------------------
612
613 *fLog << inf;
614 fLog->Separator(GetDescriptor());
615 *fLog << "Perform cuts for data set " << set.GetBaseName() << endl;
616 *fLog << endl;
617
618 // --------------------------------------------------------------------------------
619
620 // Setup Parlist
621 plist.AddToList(this); // take care of fDisplay!
622
623 MParameterI par("DataType");
624 plist.AddToList(&par);
625
626 // Setup Tasklist
627 MTaskList tlist;
628 plist.AddToList(&tlist);
629
630 // La Palma Magic1
631 MObservatory obs;
632 plist.AddToList(&obs);
633
634 if (source.IsInitialized())
635 plist.AddToList(&source);
636
637 // Initialize default binnings
638 MBinning bins1( 18, 0, 90, "BinningAlpha", "lin");
639 MBinning bins2( 15, 10, 1e6 , "BinningSize", "log");
640 MBinning bins3( 67, -0.005, 0.665, "BinningTheta", "asin");
641 //MBinning binsT(150, 0, 150, "BinningDeltaT", "lin");
642 MBinning bins4("BinningFalseSource");
643 MBinning bins5("BinningWidth");
644 MBinning bins6("BinningLength");
645 MBinning bins7("BinningDist");
646 MBinning bins8("BinningSlope");
647 MBinning bins9("BinningM3Long");
648 MBinning bins0("BinningM3Trans");
649 MBinning binsa("BinningAsym");
650 MBinning binsb("BinningConc1");
651 plist.AddToList(&bins1);
652 plist.AddToList(&bins2);
653 plist.AddToList(&bins3);
654 plist.AddToList(&bins4);
655 plist.AddToList(&bins5);
656 plist.AddToList(&bins6);
657 plist.AddToList(&bins7);
658 plist.AddToList(&bins8);
659 plist.AddToList(&bins9);
660 plist.AddToList(&bins0);
661 plist.AddToList(&binsa);
662 plist.AddToList(&binsb);
663 //plist.AddToList(&binsT);
664
665 // --------------------------------------------------------------------------------
666
667 MParameterD scale;
668 scale.SetVal(1./fNumOffSourcePos);
669
670 // Setup fitter and histograms
671 MAlphaFitter fit;
672 plist.AddToList(&fit);
673 if (set.IsWobbleMode())
674 fit.SetScaleUser(1./fNumOffSourcePos); // includes fit.SetScaleMode(MAlphaFitter::kUserScale);
675
676 MHAlpha *halphaoff = CreateNewHist(plist, "Off");
677 MFillH falpha(halphaoff, "", "FillHist");
678 MH *hfsoff = CreateNewHistFS(plist, "Off");
679 MFillH ffs(hfsoff, "MHillas", "FillFS");
680
681 // FIXME: If fPathIn read cuts and energy estimator from file!
682 MContinue contq("", "CutQ");
683 MContinue cont0("", "Cut0");
684 MContinue cont1("", "Cut1");
685 MContinue cont2("", "Cut2");
686 MContinue cont3("", "Cut3");
687 contq.SetAllowEmpty();
688 cont0.SetAllowEmpty();
689 cont1.SetAllowEmpty();
690 cont2.SetAllowEmpty();
691 cont3.SetAllowEmpty();
692
693 // Filter for VsSize
694 MFDataPhrase ftheta(0, "CutT");
695
696 // ------------- Loop Off Data --------------------
697 MReadReports readoffdata;
698 readoffdata.AddTree("Events", "MTime.", MReadReports::kMaster);
699 readoffdata.AddTree("Drive", MReadReports::kRequired);
700 readoffdata.AddTree("Starguider", MReadReports::kRequired);
701 readoffdata.AddTree("EffectiveOnTime");
702
703 MReadMarsFile readoffmc("Events");
704 readoffmc.DisableAutoScheme();
705
706 MRead &readoff = set.IsMonteCarlo() ? (MRead&)readoffmc : (MRead&)readoffdata;
707 const Bool_t setrc = set.IsWobbleMode() ? set.AddFilesOn(readoff) : set.AddFilesOff(readoff);
708 if (!setrc && (set.HasOffSequences() || set.IsWobbleMode()))
709 {
710 *fLog << err << "MDataSet::AddFiles" << (set.IsWobbleMode()?"On":"Off") << " failed." << endl;
711 return kFALSE;
712 }
713
714 const TString path(Form("%s/", fPathOut.Data()));
715 TString fname0(path);
716 TString fname1(path);
717 fname0 += fNameSummary.IsNull() ? (TString) Form("ganymed%08d-summary.root", set.GetNumAnalysis()) : fNameSummary;
718 fname1 += fNameResult.IsNull() ? (TString) Form("ganymed%08d.root", set.GetNumAnalysis()) : fNameResult;
719
720 MWriteRootFile dummy0(fPathOut.IsNull()||!CanStoreSummary()?0:fname0.Data(), fOverwrite?"RECREATE":"NEW");
721 MWriteRootFile dummy1(fPathOut.IsNull()||!CanStoreResult() ?0:fname1.Data(), fOverwrite?"RECREATE":"NEW");
722
723 MWriteRootFile *write0 = CanStoreSummary() ? &dummy0 : 0;
724 MWriteRootFile *write1 = CanStoreResult() ? &dummy1 : 0;
725 SetupWriter(write0, "WriteAfterCut0");
726 SetupWriter(write1, "WriteAfterCut3");
727
728 MTaskEnv taskenv0("CalcDisp");
729 taskenv0.SetDefault(fCalcDisp);
730
731 MTaskEnv taskenv1("CalcHadronness");
732 taskenv1.SetDefault(fCalcHadronness);
733
734 MTaskEnv taskenv2("EstimateEnergy");
735 taskenv2.SetDefault(fEstimateEnergy);
736
737 MParameterCalc setevtnum("MRawEvtHeader.fDAQEvtNumber", "SetEvtNumber");
738 setevtnum.SetNameParameter("EvtNumber");
739
740 MParameterCalc setrunnum("MRawRunHeader.GetFileID", "SetFileId");
741 setrunnum.SetNameParameter("FileId");
742
743 MFillH fill1a("MHHillasOffPre [MHHillas]", "MHillas", "FillHillasPre");
744 MFillH fill2a("MHHillasOffPost [MHHillas]", "MHillas", "FillHillasPost");
745 MFillH fill3a("MHVsSizeOffPost [MHVsSize]", "MHillasSrc", "FillVsSizePost");
746 MFillH fill3c("MHVsSizeOffTheta [MHVsSize]", "MHillasSrc", "FillVsSizeTheta");
747 MFillH fill4a("MHHilExtOffPost [MHHillasExt]", "MHillasSrc", "FillHilExtPost");
748 MFillH fill5a("MHHilSrcOffPost [MHHillasSrc]", "MHillasSrc", "FillHilSrcPost");
749 MFillH fill6a("MHImgParOffPost [MHImagePar]", "MImagePar", "FillImgParPost");
750 MFillH fill7a("MHNewParOffPost [MHNewImagePar]", "MNewImagePar", "FillNewParPost");
751 //MFillH fill9a("MHEffOffTime [MHEffectiveOnTime]", "MTime", "FillEffOnTime");
752 fill1a.SetNameTab("PreCut");
753 fill2a.SetNameTab("PostCut");
754 fill3a.SetNameTab("VsSize");
755 fill3c.SetNameTab("CutT");
756 fill4a.SetNameTab("HilExt");
757 fill5a.SetNameTab("HilSrc");
758 fill6a.SetNameTab("ImgPar");
759 fill7a.SetNameTab("NewPar");
760 //fill9a.SetNameTab("EffOffT");
761
762 fill3c.SetFilter(&ftheta);
763
764 //MFDataMember fbin("Bin.fVal", '>', 0);
765 //fill9a.SetFilter(&fbin);
766
767 MPrint print2("MEffectiveOnTime");
768 print2.EnableSkip();
769
770 MTaskList tlist2;
771 if (set.IsWobbleMode())
772 {
773 tlist2.SetNumPasses(fNumOffSourcePos);
774 fill2a.SetWeight(&scale);
775 fill3a.SetWeight(&scale);
776 fill3c.SetWeight(&scale);
777 fill4a.SetWeight(&scale);
778 fill5a.SetWeight(&scale);
779 fill6a.SetWeight(&scale);
780 fill7a.SetWeight(&scale);
781 }
782
783 // How to get source position from off- and on-data?
784 MSrcPosCalc scalc;
785 scalc.SetMode(set.IsWobbleMode()?MSrcPosCalc::kWobble:MSrcPosCalc::kOffData); /********************/
786 scalc.SetCallback(&tlist2);
787
788 MSrcPosCorrect scor;
789
790 MHillasCalc hcalc;
791 MHillasCalc hcalc2("MHillasCalcAnti");
792 hcalc.SetFlags(MHillasCalc::kCalcHillasSrc);
793 hcalc2.SetFlags(MHillasCalc::kCalcHillasSrc);
794 hcalc2.SetNameHillasSrc("MHillasSrcAnti");
795 hcalc2.SetNameSrcPosCam("MSrcPosAnti");
796
797 MSrcPosRndm srcrndm;
798
799 MH3 hvs("MPointingPos.fZd");
800 hvs.SetName("ThetaOff;Theta");
801 hvs.SetTitle("Effective On-Time vs. Zenith Angle;\\Theta [\\circ];T_{on} [s]");
802
803 MFillH fillvs(&hvs, "", "FillOnTime");
804 if (!set.IsMonteCarlo())
805 fillvs.SetWeight("MEffectiveOnTime");
806 fillvs.SetNameTab("OnTime");
807
808 // It is not really necessary to re-calculate the image parameters
809 // for the the on-source for MCs, but it is done for symmetry reasons
810 if (set.IsMonteCarlo())
811 tlist2.AddToList(&fillvs);
812
813 tlist2.AddToList(&scalc);
814 tlist2.AddToList(&scor);
815 tlist2.AddToList(&srcrndm);
816 tlist2.AddToList(&hcalc);
817 if (set.IsWobbleMode())
818 tlist2.AddToList(&hcalc2);
819 tlist2.AddToList(&cont0);
820 tlist2.AddToList(&taskenv0);
821 tlist2.AddToList(&taskenv1);
822 tlist2.AddToList(&setrunnum);
823 tlist2.AddToList(&setevtnum);
824 if (write0)
825 tlist2.AddToList(write0);
826 tlist2.AddToList(&cont1);
827 if (!fWriteOnly && (!set.IsWobbleMode() || !fNameHistFS.IsNull()))
828 tlist2.AddToList(&ffs);
829 tlist2.AddToList(&cont2);
830 if (!fWriteOnly)
831 {
832 tlist2.AddToList(&fill2a);
833 if (fFullDisplay)
834 {
835 tlist2.AddToList(&ftheta);
836 tlist2.AddToList(&fill3a);
837 tlist2.AddToList(&fill3c);
838 tlist2.AddToList(&fill4a);
839 tlist2.AddToList(&fill5a);
840 tlist2.AddToList(&fill6a);
841 tlist2.AddToList(&fill7a);
842 }
843 }
844 if (!fWriteOnly)
845 {
846 tlist2.AddToList(&falpha);
847 /* if (!fIsMonteCarlo)
848 {
849 tlist2.AddToList(&fbin);
850 tlist2.AddToList(&fill9a);
851 } */
852 }
853 tlist2.AddToList(&cont3);
854 tlist2.AddToList(&taskenv2);
855
856 if (write1)
857 tlist2.AddToList(write1);
858
859 MPointingDevCalc devcalc;
860
861 tlist.AddToList(&readoff);
862 if (!set.IsMonteCarlo())
863 tlist.AddToList(&fillvs, "EffectiveOnTime");
864 if (gLog.GetDebugLevel()>4)
865 tlist.AddToList(&print2, "EffectiveOnTime");
866 tlist.AddToList(&devcalc, "Starguider");
867 tlist.AddToList(&contq, "Events");
868 if (!fWriteOnly)
869 tlist.AddToList(&fill1a, "Events");
870 tlist.AddToList(&tlist2, "Events");
871
872 // by setting it here it is distributed to all consecutive tasks
873 tlist.SetAccelerator(MTask::kAccDontReset|MTask::kAccDontTime);
874
875 par.SetVal(0);
876
877 // Create and setup the eventloop
878 MEvtLoop evtloop(fName);
879 evtloop.SetParList(&plist);
880 evtloop.SetDisplay(fDisplay);
881 evtloop.SetLogStream(fLog);
882 if (!SetupEnv(evtloop))
883 return kFALSE;
884
885 TObjArray cont;
886 cont.Add(&contq);
887 cont.Add(&cont0);
888 cont.Add(&cont1);
889 cont.Add(&cont2);
890 cont.Add(&cont3);
891 //if (taskenv1.GetTask())
892 // cont.Add(taskenv1.GetTask());
893 if (taskenv0.GetTask())
894 cont.Add(taskenv0.GetTask());
895 if (taskenv1.GetTask())
896 cont.Add(taskenv1.GetTask());
897 if (taskenv2.GetTask())
898 cont.Add(taskenv2.GetTask());
899
900 if (!WriteTasks(set.GetNumAnalysis(), cont))
901 return kFALSE;
902
903 if (set.HasOffSequences() || set.IsWobbleMode())
904 {
905 // Execute first analysis
906 if (!evtloop.Eventloop(fMaxEvents))
907 {
908 *fLog << err << GetDescriptor() << ": Processing of off-sequences failed." << endl;
909 return -3;
910 }
911
912 if (!evtloop.GetDisplay())
913 {
914 *fLog << err << GetDescriptor() << ": Execution stopped by user." << endl;
915 return kFALSE;
916 }
917
918 //plist.FindObject("MTimeEffectiveOnTime")->Clear();
919 }
920 else
921 {
922 // This is the simplest way to remove the two object from the parlist
923 delete halphaoff;
924 delete hfsoff;
925 }
926
927 // ------------- Loop On Data --------------------
928 MReadReports readondata;
929 readondata.AddTree("Events", "MTime.", MReadReports::kMaster);
930 readondata.AddTree("Drive", MReadReports::kRequired);
931 readondata.AddTree("Starguider", MReadReports::kRequired);
932 readondata.AddTree("EffectiveOnTime");
933
934 MReadMarsFile readonmc("Events");
935 readonmc.DisableAutoScheme();
936
937 MRead &readon = set.IsMonteCarlo() ? (MRead&)readonmc : (MRead&)readondata;
938 if (!set.AddFilesOn(readon))
939 return kFALSE;
940
941 scalc.SetMode(MSrcPosCalc::kDefault);
942 scalc.SetNumRandomOffPositions(fNumOffSourcePos);
943
944 MFillH fill1b("MHHillasOnPre [MHHillas]", "MHillas", "FillHillasPre");
945 MFillH fill2b("MHHillasOnPost [MHHillas]", "MHillas", "FillHillasPost");
946 MFillH fill3b("MHVsSizeOnPost [MHVsSize]", "MHillasSrc", "FillVsSizePost");
947 MFillH fill3d("MHVsSizeOnTheta [MHVsSize]", "MHillasSrc", "FillVsSizeTheta");
948 MFillH fill4b("MHHilExtOnPost [MHHillasExt]", "MHillasSrc", "FillHilExtPost");
949 MFillH fill5b("MHHilSrcOnPost [MHHillasSrc]", "MHillasSrc", "FillHilSrcPost");
950 MFillH fill6b("MHImgParOnPost [MHImagePar]", "MImagePar", "FillImgParPost");
951 MFillH fill7b("MHNewParOnPost [MHNewImagePar]", "MNewImagePar", "FillNewParPost");
952 //MFillH fill9b("MHEffOnTime [MHEffectiveOnTime]", "MTime", "FillEffOnTime");
953 fill1b.SetNameTab("PreCut");
954 fill2b.SetNameTab("PostCut");
955 fill3b.SetNameTab("VsSize");
956 fill3d.SetNameTab("CutT");
957 fill4b.SetNameTab("HilExt");
958 fill5b.SetNameTab("HilSrc");
959 fill6b.SetNameTab("ImgPar");
960 fill7b.SetNameTab("NewPar");
961 //fill9b.SetNameTab("EffOnT");
962 fill1b.SetDrawOption(set.HasOffSequences()||set.IsWobbleMode()?"same":"");
963 fill2b.SetDrawOption(set.HasOffSequences()||set.IsWobbleMode()?"same":"");
964 fill3b.SetDrawOption(set.HasOffSequences()||set.IsWobbleMode()?"same":"");
965 fill3d.SetDrawOption(set.HasOffSequences()||set.IsWobbleMode()?"same":"");
966 fill4b.SetDrawOption(set.HasOffSequences()||set.IsWobbleMode()?"same":"");
967 fill5b.SetDrawOption(set.HasOffSequences()||set.IsWobbleMode()?"same":"");
968 fill6b.SetDrawOption(set.HasOffSequences()||set.IsWobbleMode()?"same":"");
969 fill7b.SetDrawOption(set.HasOffSequences()||set.IsWobbleMode()?"same":"");
970 //fill9b.SetFilter(&fbin);
971
972 fill3d.SetFilter(&ftheta);
973
974 /*
975 MHVsTime hvs("MEffectiveOnTime.fVal");
976 hvs.SetTitle("Effective On-Time vs. Time;;T_{on}");
977 MFillH fillvs(&hvs, "MTimeEffectiveOnTime", "FillOnTime");
978 fillvs.SetNameTab("OnTime");
979 */
980
981 /*
982 MParameterD weight;
983 weight.SetVal(-1);
984 fill2a.SetWeight(&weight);
985 fill3a.SetWeight(&weight);
986 fill4a.SetWeight(&weight);
987 fill5a.SetWeight(&weight);
988 fill6a.SetWeight(&weight);
989 fill7a.SetWeight(&weight);
990 if (fSubstraction)
991 {
992 fill2a.SetNameTab("PostCut-");
993 fill3a.SetNameTab("VsSize-");
994 fill4a.SetNameTab("HilExt-");
995 fill5a.SetNameTab("HilSrc-");
996 fill6a.SetNameTab("ImgPar-");
997 fill7a.SetNameTab("NewPar-");
998 }
999 */
1000 MHAlpha *halphaon=CreateNewHist(plist);
1001 MFillH falpha2(halphaon, "", "FillHist");
1002 MH *hfs=CreateNewHistFS(plist);
1003 MFillH ffs2(hfs, "MHillas", "FillFS");
1004 MFillH fillphi("MHPhi", "", "FillPhi");
1005 fillphi.SetDrawOption("anticut");
1006
1007 tlist2.SetNumPasses();
1008
1009 tlist.Replace(&readon);
1010 if (fRndmSrcPos && !set.IsWobbleMode())
1011 tlist2.RemoveFromList(&srcrndm);
1012
1013 MFillH fillsrc(&hsrcpos, "MSrcPosCam", "FillSrcPosCam");
1014 fillsrc.SetNameTab("SrcPos");
1015
1016 if (set.IsWobbleMode()/* && !set.IsMonteCarlo()*/)
1017 tlist2.AddToListBefore(&fillsrc, &hcalc);
1018
1019 MH3 hvs2("MPointingPos.fZd");
1020 hvs2.SetName("Theta");
1021 hvs2.SetTitle("Effective On-Time vs. Zenith Angle;\\Theta [\\circ];T_{on} [s]");
1022
1023 MFillH fillvs2(&hvs2, "", "FillOnTime");
1024 if (!set.IsMonteCarlo())
1025 fillvs2.SetWeight("MEffectiveOnTime");
1026 fillvs2.SetNameTab("OnTime");
1027 fillvs2.SetDrawOption(set.HasOffSequences()||set.IsWobbleMode()?"same":"");
1028
1029 if (!fWriteOnly)
1030 {
1031 tlist.Replace(&fill1b);
1032
1033 tlist2.Replace(&fill2b);
1034 if (fFullDisplay)
1035 {
1036 tlist2.Replace(&fill3b);
1037 tlist2.Replace(&fill3d);
1038 tlist2.Replace(&fill4b);
1039 tlist2.Replace(&fill5b);
1040 tlist2.Replace(&fill6b);
1041 tlist2.Replace(&fill7b);
1042 }
1043 tlist2.Replace(&falpha2);
1044 //if (!fIsMonteCarlo)
1045 // tlist2.Replace(&fill9b);
1046 if (!set.IsWobbleMode() || !fNameHist.IsNull())
1047 tlist2.Replace(&ffs2);
1048 if (set.IsWobbleMode())
1049 {
1050 tlist2.AddToListAfter(&fillphi, &falpha2);
1051 if (!fNameHist.IsNull())
1052 tlist2.RemoveFromList(&ffs);
1053 }
1054
1055 if (!set.IsMonteCarlo())
1056 tlist.Replace(&fillvs2);
1057 else
1058 tlist2.Replace(&fillvs2);
1059 }
1060
1061 // by setting it here it is distributed to all consecutive tasks
1062 tlist.SetAccelerator(MTask::kAccDontReset|MTask::kAccDontTime);
1063
1064 par.SetVal(1);
1065
1066 // Execute first analysis
1067 if (!evtloop.Eventloop(fMaxEvents))
1068 {
1069 *fLog << err << GetDescriptor() << ": Processing of on-sequences failed." << endl;
1070 return -4;
1071 }
1072
1073 // FIXME: Perform fit and plot energy dependant alpha plots
1074 // and fit result to new tabs!
1075 if (!WriteResult(plist, set))
1076 return kFALSE;
1077
1078 *fLog << all << GetDescriptor() << ": Done." << endl;
1079 *fLog << endl << endl;
1080
1081 return kTRUE;
1082}
Note: See TracBrowser for help on using the repository browser.