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

Last change on this file since 9434 was 9374, checked in by tbretz, 16 years ago
*** empty log message ***
File size: 33.9 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 const Int_t num = set.GetNumAnalysis();
324 const TString fname(fNameOutput.IsNull() ? Form("ganymed%08d.root", num) : fNameOutput.Data());
325
326 // If requested, write to already open output file
327 if (fNameResult.IsNull() && fStoreResult)
328 {
329 TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject(fname);
330 if (file)
331 {
332 file->cd();
333 return WriteContainer(arr);
334 }
335 }
336
337 return WriteContainer(arr, fname, "UPDATE");
338}
339
340// --------------------------------------------------------------------------
341//
342// MJCut allows to setup several option by a resource file:
343// MJCut.WriteSummary: yes, no
344// MJCut.SummaryFile: filename
345// MJCut.WriteResult: yes, no
346// MJCut.ResultFile: filename
347// MJCut.HistName: MHAlpha
348//
349Bool_t MJCut::CheckEnvLocal()
350{
351 const TString f0(GetEnv("SummaryFile", ""));
352 const TString f1(GetEnv("ResultFile", ""));
353 if (!f0.IsNull())
354 SetNameSummaryFile(f0);
355 if (!f1.IsNull())
356 SetNameResultFile(f1);
357
358 EnableStorageOfSummary(GetEnv("SummaryFile", fStoreSummary));
359 EnableStorageOfResult(GetEnv("ResultFile", fStoreResult));
360 EnableFullDisplay(GetEnv("FullDisplay", fFullDisplay));
361 EnableRandomSrcPos(GetEnv("RandomSourcePosition", fRndmSrcPos));
362
363 fNumOffSourcePos = GetEnv("NumOffSourcePositions", (Int_t)fNumOffSourcePos);
364 //EnableSubstraction(GetEnv("HistogramSubstraction", fSubstraction));
365
366 SetNameHist(GetEnv("NameHist", fNameHist));
367 SetNameHistFS(GetEnv("NameHistFS", fNameHistFS));
368
369 return kTRUE;
370}
371
372// --------------------------------------------------------------------------
373//
374// Setup write to write:
375// container tree optional?
376// -------------- ---------- -----------
377// "MHillas" to "Events"
378// "MHillasSrc" to "Events"
379// "Hadronness" to "Events" yes
380// "MEnergyEst" to "Events" yes
381// "DataType" to "Events"
382//
383void MJCut::SetupWriter(MWriteRootFile *write, const char *name) const
384{
385 if (!write)
386 return;
387
388 write->SetName(name);
389 write->AddContainer("MHillas", "Events");
390 write->AddContainer("MHillasSrc", "Events");
391 write->AddContainer("MHillasExt", "Events");
392 write->AddContainer("MPointingPos", "Events");
393 write->AddContainer("MHillasSrcAnti", "Events", kFALSE);
394 write->AddContainer("MImagePar", "Events", kFALSE);
395 write->AddContainer("MNewImagePar", "Events", kFALSE);
396 write->AddContainer("MNewImagePar2", "Events", kFALSE);
397 write->AddContainer("Hadronness", "Events", kFALSE);
398 write->AddContainer("MSrcPosCam", "Events", kFALSE);
399 write->AddContainer("MSrcPosAnti", "Events", kFALSE);
400 write->AddContainer("ThetaSquared", "Events", kFALSE);
401 write->AddContainer("OpticalAxis", "Events", kFALSE);
402 write->AddContainer("Disp", "Events", kFALSE);
403 write->AddContainer("Ghostbuster", "Events", kFALSE);
404 write->AddContainer("MEnergyEst", "Events", kFALSE);
405 write->AddContainer("MTime", "Events", kFALSE);
406 write->AddContainer("MMcEvt", "Events", kFALSE);
407 write->AddContainer("DataType", "Events");
408 write->AddContainer("FileId", "Events");
409 write->AddContainer("EvtNumber", "Events");
410 // write->AddContainer("MMuonSearchPar", "Events", kFALSE);
411 // write->AddContainer("MMuonCalibPar", "Events", kFALSE);
412}
413
414// --------------------------------------------------------------------------
415//
416// Create a new instance of an object with name name of class
417// type fNameHist in parlist. It must derive from MHAlpha.
418// Call ForceUsingSize for it and return its pointer.
419// If something fails NULL is returned.
420//
421MHAlpha *MJCut::CreateNewHist(MParList &plist, const char *name) const
422{
423 TClass *cls = gROOT->GetClass(fNameHist);
424 if (!cls)
425 {
426 *fLog << err << "Class " << fNameHist << " not found in dictionary... abort." << endl;
427 return NULL;
428 }
429 if (!cls->InheritsFrom(MHAlpha::Class()))
430 {
431 *fLog << err << "Class " << fNameHist << " doesn't inherit from MHAlpha... abort." << endl;
432 return NULL;
433 }
434
435 const TString objname(Form("Hist%s", name));
436 MHAlpha *h = (MHAlpha*)plist.FindCreateObj(fNameHist, objname);
437 if (!h)
438 return NULL;
439
440 h->ForceUsingSize();
441
442 return h;
443}
444
445// --------------------------------------------------------------------------
446//
447// Create a new instance of an object with name name of class
448// type fNameHistFS in parlist. It must derive from MHFalseSource
449// If something fails NULL is returned.
450//
451MH *MJCut::CreateNewHistFS(MParList &plist, const char *name) const
452{
453 const TString cname(fNameHistFS.IsNull()?"MHFalseSource":fNameHistFS.Data());
454
455 TClass *cls = gROOT->GetClass(cname);
456 if (!cls)
457 {
458 *fLog << err << "Class " << cname << " not found in dictionary... abort." << endl;
459 return NULL;
460 }
461 if (!cls->InheritsFrom("MHFalseSource"))
462 {
463 *fLog << err << "Class " << cname << " doesn't inherit from MHFalseSource... abort." << endl;
464 return NULL;
465 }
466
467 const TString objname(Form("FS%s", name));
468 return (MH*)plist.FindCreateObj(cname, objname);
469}
470
471Bool_t MJCut::FillSrcPosCam(const MDataSet &set, MPointingPos &source, MHSrcPosCam &hsrcpos)
472{
473 *fLog << inf;
474 fLog->Separator(GetDescriptor());
475 *fLog << "Filling MHSrcPosCam " << set.GetBaseName() << endl;
476 *fLog << endl;
477
478 // --------------------------------------------------------------------------------
479
480 // Setup Parlist
481 MParList plist;
482 plist.AddToList(this); // take care of fDisplay!
483
484 // Setup Tasklist
485 MTaskList tlist;
486 plist.AddToList(&tlist);
487
488 // La Palma Magic1, Possible source position
489 MObservatory obs;
490 plist.AddToList(&obs);
491 plist.AddToList(&source);
492
493 // Initialize default binnings
494 // MBinning bins1(18, 0, 90, "BinningSrcPosCam", "lin");
495 // plist.AddToList(&bins1);
496
497 // ------------- Loop Off Data --------------------
498 MReadReports read;
499
500 read.EnableAutoScheme();
501 read.AddToBranchList("MTimeEffectiveOnTime.*");
502 read.AddToBranchList("MEffectiveOnTime.*");
503
504 read.AddTree("Events", "MTime.", MReadReports::kMaster);
505 read.AddTree("Drive", MReadReports::kRequired);
506 read.AddTree("Starguider", MReadReports::kRequired);
507 read.AddTree("EffectiveOnTime");
508
509 if (!set.AddFilesOn(read))
510 return kFALSE;
511
512 MFillH fill(&hsrcpos, "MSrcPosCam", "FillSrcPosCam");
513 fill.SetNameTab("SrcPos");
514
515 // How to get source position from off- and on-data?
516 MSrcPosCorrect scor;
517 MSrcPosCalc scalc;
518 scalc.SetMode(MSrcPosCalc::kDefault);
519
520 MPointingDevCalc devcalc;
521
522 tlist.AddToList(&read);
523 tlist.AddToList(&devcalc, "Starguider");
524 tlist.AddToList(&scalc, "Events");
525 tlist.AddToList(&scor, "Events");
526 tlist.AddToList(&fill, "Events");
527
528 // by setting it here it is distributed to all consecutive tasks
529 tlist.SetAccelerator(MTask::kAccDontReset|MTask::kAccDontTime);
530
531 // Create and setup the eventloop
532 MEvtLoop evtloop(fName);
533 evtloop.SetParList(&plist);
534 evtloop.SetDisplay(fDisplay);
535 evtloop.SetLogStream(fLog);
536 if (!SetupEnv(evtloop))
537 return kFALSE;
538
539 // Execute first analysis
540 if (!evtloop.Eventloop(fMaxEvents))
541 {
542 *fLog << err << GetDescriptor() << ": Processing of on-sequences to fill SrcPosCam failed." << endl;
543 return kFALSE;
544 }
545
546 if (!evtloop.GetDisplay())
547 {
548 *fLog << err << GetDescriptor() << ": Execution stopped by user." << endl;
549 return kFALSE;
550 }
551
552 *fLog << all << GetDescriptor() << ": Done." << endl;
553 *fLog << endl << endl;
554
555 return kTRUE;
556}
557
558Int_t MJCut::Process(const MDataSet &set)
559{
560 if (!set.IsValid())
561 {
562 *fLog << err << "ERROR - DataSet invalid!" << endl;
563 return kFALSE;
564 }
565
566 if (!HasWritePermission(GetOutputFile(set.GetNumAnalysis())))
567 return kFALSE;
568
569 if (!CheckEnv())
570 return kFALSE;
571
572 // --------------------------------------------------------------------------------
573
574 if (!set.IsWobbleMode() && fNumOffSourcePos!=1)
575 {
576 *fLog << inf << "No wobble mode, but NumOffSoucePos!=1 (" << fNumOffSourcePos << ")... reset to 1." << endl;
577 fNumOffSourcePos = 1;
578 }
579
580 // Possible source position (eg. Wobble Mode)
581 MPointingPos source("MSourcePos");
582 if (set.HasSource())
583 {
584 if (!set.GetSourcePos(source))
585 return -1;
586 *fLog << all;
587 source.Print("RaDec");
588 }
589 else
590 *fLog << all << "No source position applied..." << endl;
591
592 // ADD A CHECK FOR THE SOURCE POSITION by COMPARISON WITH RA/DEC
593
594 MParList plist;
595
596 MHSrcPosCam hsrcpos(set.IsWobbleMode());
597 if (!set.IsWobbleMode() && source.IsInitialized() && fRndmSrcPos)
598 {
599 if (!FillSrcPosCam(set, source, hsrcpos))
600 return -2;
601 plist.AddToList(&hsrcpos);
602 }
603
604 // --------------------------------------------------------------------------------
605
606 *fLog << inf;
607 fLog->Separator(GetDescriptor());
608 *fLog << "Perform cuts for data set " << set.GetBaseName() << endl;
609 *fLog << endl;
610
611 // --------------------------------------------------------------------------------
612
613 // Setup Parlist
614 plist.AddToList(this); // take care of fDisplay!
615
616 MParameterI par("DataType");
617 plist.AddToList(&par);
618
619 // Setup Tasklist
620 MTaskList tlist;
621 plist.AddToList(&tlist);
622
623 // La Palma Magic1
624 MObservatory obs;
625 plist.AddToList(&obs);
626
627 if (source.IsInitialized())
628 plist.AddToList(&source);
629
630 // Initialize default binnings
631 MBinning bins1( 18, 0, 90, "BinningAlpha", "lin");
632 MBinning bins2( 15, 10, 1e6 , "BinningSize", "log");
633 MBinning bins3( 67, -0.005, 0.665, "BinningTheta", "asin");
634 //MBinning binsT(150, 0, 150, "BinningDeltaT", "lin");
635 MBinning bins4("BinningFalseSource");
636 MBinning bins5("BinningWidth");
637 MBinning bins6("BinningLength");
638 MBinning bins7("BinningDist");
639 MBinning bins8("BinningSlope");
640 MBinning bins9("BinningM3Long");
641 MBinning bins0("BinningM3Trans");
642 MBinning binsa("BinningAsym");
643 MBinning binsb("BinningConc1");
644 plist.AddToList(&bins1);
645 plist.AddToList(&bins2);
646 plist.AddToList(&bins3);
647 plist.AddToList(&bins4);
648 plist.AddToList(&bins5);
649 plist.AddToList(&bins6);
650 plist.AddToList(&bins7);
651 plist.AddToList(&bins8);
652 plist.AddToList(&bins9);
653 plist.AddToList(&bins0);
654 plist.AddToList(&binsa);
655 plist.AddToList(&binsb);
656 //plist.AddToList(&binsT);
657
658 // --------------------------------------------------------------------------------
659
660 MParameterD scale;
661 scale.SetVal(1./fNumOffSourcePos);
662
663 // Setup fitter and histograms
664 MAlphaFitter fit;
665 plist.AddToList(&fit);
666 if (set.IsWobbleMode())
667 fit.SetScaleUser(1./fNumOffSourcePos); // includes fit.SetScaleMode(MAlphaFitter::kUserScale);
668
669 MHAlpha *halphaoff = CreateNewHist(plist, "Off");
670 MFillH falpha(halphaoff, "", "FillHist");
671 MH *hfsoff = CreateNewHistFS(plist, "Off");
672 MFillH ffs(hfsoff, "MHillas", "FillFS");
673
674 // FIXME: If fPathIn read cuts and energy estimator from file!
675 MContinue contq("", "CutQ");
676 MContinue cont0("", "Cut0");
677 MContinue cont1("", "Cut1");
678 MContinue cont2("", "Cut2");
679 MContinue cont3("", "Cut3");
680 contq.SetAllowEmpty();
681 cont0.SetAllowEmpty();
682 cont1.SetAllowEmpty();
683 cont2.SetAllowEmpty();
684 cont3.SetAllowEmpty();
685
686 // Filter for VsSize
687 MFDataPhrase ftheta(0, "CutT");
688
689 // ------------- Loop Off Data --------------------
690 MReadReports readoffdata;
691 readoffdata.AddTree("Events", "MTime.", MReadReports::kMaster);
692 readoffdata.AddTree("Drive", MReadReports::kRequired);
693 readoffdata.AddTree("Starguider", MReadReports::kRequired);
694 readoffdata.AddTree("EffectiveOnTime");
695
696 MReadMarsFile readoffmc("Events");
697 readoffmc.DisableAutoScheme();
698
699 MRead &readoff = set.IsMonteCarlo() ? (MRead&)readoffmc : (MRead&)readoffdata;
700 const Bool_t setrc = set.IsWobbleMode() ? set.AddFilesOn(readoff) : set.AddFilesOff(readoff);
701 if (!setrc && (set.HasOffSequences() || set.IsWobbleMode()))
702 {
703 *fLog << err << "MDataSet::AddFiles" << (set.IsWobbleMode()?"On":"Off") << " failed." << endl;
704 return kFALSE;
705 }
706
707 const TString path(Form("%s/", fPathOut.Data()));
708 TString fname0(path);
709 TString fname1(path);
710 fname0 += fNameSummary.IsNull() ? (TString) Form("ganymed%08d-summary.root", set.GetNumAnalysis()) : fNameSummary;
711 fname1 += fNameResult.IsNull() ? (TString) Form("ganymed%08d.root", set.GetNumAnalysis()) : fNameResult;
712
713 MWriteRootFile dummy0(fPathOut.IsNull()||!CanStoreSummary()?0:fname0.Data(), fOverwrite?"RECREATE":"NEW");
714 MWriteRootFile dummy1(fPathOut.IsNull()||!CanStoreResult() ?0:fname1.Data(), fOverwrite?"RECREATE":"NEW");
715
716 MWriteRootFile *write0 = CanStoreSummary() ? &dummy0 : 0;
717 MWriteRootFile *write1 = CanStoreResult() ? &dummy1 : 0;
718 SetupWriter(write0, "WriteAfterCut0");
719 SetupWriter(write1, "WriteAfterCut3");
720
721 MTaskEnv taskenv0("CalcDisp");
722 taskenv0.SetDefault(fCalcDisp);
723
724 MTaskEnv taskenv1("CalcHadronness");
725 taskenv1.SetDefault(fCalcHadronness);
726
727 MTaskEnv taskenv2("EstimateEnergy");
728 taskenv2.SetDefault(fEstimateEnergy);
729
730 MParameterCalc setevtnum("MRawEvtHeader.fDAQEvtNumber", "SetEvtNumber");
731 setevtnum.SetNameParameter("EvtNumber");
732
733 MParameterCalc setrunnum("MRawRunHeader.GetFileID", "SetFileId");
734 setrunnum.SetNameParameter("FileId");
735
736 MFillH fill1a("MHHillasOffPre [MHHillas]", "MHillas", "FillHillasPre");
737 MFillH fill2a("MHHillasOffPost [MHHillas]", "MHillas", "FillHillasPost");
738 MFillH fill3a("MHVsSizeOffPost [MHVsSize]", "MHillasSrc", "FillVsSizePost");
739 MFillH fill3c("MHVsSizeOffTheta [MHVsSize]", "MHillasSrc", "FillVsSizeTheta");
740 MFillH fill4a("MHHilExtOffPost [MHHillasExt]", "MHillasSrc", "FillHilExtPost");
741 MFillH fill5a("MHHilSrcOffPost [MHHillasSrc]", "MHillasSrc", "FillHilSrcPost");
742 MFillH fill6a("MHImgParOffPost [MHImagePar]", "MImagePar", "FillImgParPost");
743 MFillH fill7a("MHNewParOffPost [MHNewImagePar]", "MNewImagePar", "FillNewParPost");
744 //MFillH fill9a("MHEffOffTime [MHEffectiveOnTime]", "MTime", "FillEffOnTime");
745 fill1a.SetNameTab("PreCut");
746 fill2a.SetNameTab("PostCut");
747 fill3a.SetNameTab("VsSize");
748 fill3c.SetNameTab("CutT");
749 fill4a.SetNameTab("HilExt");
750 fill5a.SetNameTab("HilSrc");
751 fill6a.SetNameTab("ImgPar");
752 fill7a.SetNameTab("NewPar");
753 //fill9a.SetNameTab("EffOffT");
754
755 fill3c.SetFilter(&ftheta);
756
757 //MFDataMember fbin("Bin.fVal", '>', 0);
758 //fill9a.SetFilter(&fbin);
759
760 MPrint print2("MEffectiveOnTime");
761 print2.EnableSkip();
762
763 MTaskList tlist2;
764 if (set.IsWobbleMode())
765 {
766 tlist2.SetNumPasses(fNumOffSourcePos);
767 fill2a.SetWeight(&scale);
768 fill3a.SetWeight(&scale);
769 fill3c.SetWeight(&scale);
770 fill4a.SetWeight(&scale);
771 fill5a.SetWeight(&scale);
772 fill6a.SetWeight(&scale);
773 fill7a.SetWeight(&scale);
774 }
775
776 // How to get source position from off- and on-data?
777 MSrcPosCalc scalc;
778 scalc.SetMode(set.IsWobbleMode()?MSrcPosCalc::kWobble:MSrcPosCalc::kOffData); /********************/
779 scalc.SetCallback(&tlist2);
780
781 MSrcPosCorrect scor;
782
783 MHillasCalc hcalc;
784 MHillasCalc hcalc2("MHillasCalcAnti");
785 hcalc.SetFlags(MHillasCalc::kCalcHillasSrc);
786 hcalc2.SetFlags(MHillasCalc::kCalcHillasSrc);
787 hcalc2.SetNameHillasSrc("MHillasSrcAnti");
788 hcalc2.SetNameSrcPosCam("MSrcPosAnti");
789
790 MSrcPosRndm srcrndm;
791
792 MH3 hvs("MPointingPos.fZd");
793 hvs.SetName("ThetaOff;Theta");
794 hvs.SetTitle("Effective On-Time vs. Zenith Angle;\\Theta [\\circ];T_{on} [s]");
795
796 MFillH fillvs(&hvs, "", "FillOnTime");
797 if (!set.IsMonteCarlo())
798 fillvs.SetWeight("MEffectiveOnTime");
799 fillvs.SetNameTab("OnTime");
800
801 // It is not really necessary to re-calculate the image parameters
802 // for the the on-source for MCs, but it is done for symmetry reasons
803 if (set.IsMonteCarlo())
804 tlist2.AddToList(&fillvs);
805
806 tlist2.AddToList(&scalc);
807 tlist2.AddToList(&scor);
808 tlist2.AddToList(&srcrndm);
809 tlist2.AddToList(&hcalc);
810 if (set.IsWobbleMode())
811 tlist2.AddToList(&hcalc2);
812 tlist2.AddToList(&cont0);
813 tlist2.AddToList(&taskenv0);
814 tlist2.AddToList(&taskenv1);
815 tlist2.AddToList(&setrunnum);
816 tlist2.AddToList(&setevtnum);
817 if (write0)
818 tlist2.AddToList(write0);
819 tlist2.AddToList(&cont1);
820 if (!fWriteOnly && (!set.IsWobbleMode() || !fNameHistFS.IsNull()))
821 tlist2.AddToList(&ffs);
822 tlist2.AddToList(&cont2);
823 if (!fWriteOnly)
824 {
825 tlist2.AddToList(&fill2a);
826 if (fFullDisplay)
827 {
828 tlist2.AddToList(&ftheta);
829 tlist2.AddToList(&fill3a);
830 tlist2.AddToList(&fill3c);
831 tlist2.AddToList(&fill4a);
832 tlist2.AddToList(&fill5a);
833 tlist2.AddToList(&fill6a);
834 tlist2.AddToList(&fill7a);
835 }
836 }
837 if (!fWriteOnly)
838 {
839 tlist2.AddToList(&falpha);
840 /* if (!fIsMonteCarlo)
841 {
842 tlist2.AddToList(&fbin);
843 tlist2.AddToList(&fill9a);
844 } */
845 }
846 tlist2.AddToList(&cont3);
847 tlist2.AddToList(&taskenv2);
848
849 if (write1)
850 tlist2.AddToList(write1);
851
852 MPointingDevCalc devcalc;
853
854 tlist.AddToList(&readoff);
855 if (!set.IsMonteCarlo())
856 tlist.AddToList(&fillvs, "EffectiveOnTime");
857 if (gLog.GetDebugLevel()>4)
858 tlist.AddToList(&print2, "EffectiveOnTime");
859 tlist.AddToList(&devcalc, "Starguider");
860 tlist.AddToList(&contq, "Events");
861 if (!fWriteOnly)
862 tlist.AddToList(&fill1a, "Events");
863 tlist.AddToList(&tlist2, "Events");
864
865 // by setting it here it is distributed to all consecutive tasks
866 tlist.SetAccelerator(MTask::kAccDontReset|MTask::kAccDontTime);
867
868 par.SetVal(0);
869
870 // Create and setup the eventloop
871 MEvtLoop evtloop(fName);
872 evtloop.SetParList(&plist);
873 evtloop.SetDisplay(fDisplay);
874 evtloop.SetLogStream(fLog);
875 if (!SetupEnv(evtloop))
876 return kFALSE;
877
878 TObjArray cont;
879 cont.Add(&contq);
880 cont.Add(&cont0);
881 cont.Add(&cont1);
882 cont.Add(&cont2);
883 cont.Add(&cont3);
884 //if (taskenv1.GetTask())
885 // cont.Add(taskenv1.GetTask());
886 if (taskenv0.GetTask())
887 cont.Add(taskenv0.GetTask());
888 if (taskenv1.GetTask())
889 cont.Add(taskenv1.GetTask());
890 if (taskenv2.GetTask())
891 cont.Add(taskenv2.GetTask());
892
893 if (!WriteTasks(set.GetNumAnalysis(), cont))
894 return kFALSE;
895
896 if (set.HasOffSequences() || set.IsWobbleMode())
897 {
898 // Execute first analysis
899 if (!evtloop.Eventloop(fMaxEvents))
900 {
901 *fLog << err << GetDescriptor() << ": Processing of off-sequences failed." << endl;
902 return -3;
903 }
904
905 if (!evtloop.GetDisplay())
906 {
907 *fLog << err << GetDescriptor() << ": Execution stopped by user." << endl;
908 return kFALSE;
909 }
910
911 //plist.FindObject("MTimeEffectiveOnTime")->Clear();
912 }
913 else
914 {
915 // This is the simplest way to remove the two object from the parlist
916 delete halphaoff;
917 delete hfsoff;
918 }
919
920 // ------------- Loop On Data --------------------
921 MReadReports readondata;
922 readondata.AddTree("Events", "MTime.", MReadReports::kMaster);
923 readondata.AddTree("Drive", MReadReports::kRequired);
924 readondata.AddTree("Starguider", MReadReports::kRequired);
925 readondata.AddTree("EffectiveOnTime");
926
927 MReadMarsFile readonmc("Events");
928 readonmc.DisableAutoScheme();
929
930 MRead &readon = set.IsMonteCarlo() ? (MRead&)readonmc : (MRead&)readondata;
931 if (!set.AddFilesOn(readon))
932 return kFALSE;
933
934 scalc.SetMode(MSrcPosCalc::kDefault);
935 scalc.SetNumRandomOffPositions(fNumOffSourcePos);
936
937 MFillH fill1b("MHHillasOnPre [MHHillas]", "MHillas", "FillHillasPre");
938 MFillH fill2b("MHHillasOnPost [MHHillas]", "MHillas", "FillHillasPost");
939 MFillH fill3b("MHVsSizeOnPost [MHVsSize]", "MHillasSrc", "FillVsSizePost");
940 MFillH fill3d("MHVsSizeOnTheta [MHVsSize]", "MHillasSrc", "FillVsSizeTheta");
941 MFillH fill4b("MHHilExtOnPost [MHHillasExt]", "MHillasSrc", "FillHilExtPost");
942 MFillH fill5b("MHHilSrcOnPost [MHHillasSrc]", "MHillasSrc", "FillHilSrcPost");
943 MFillH fill6b("MHImgParOnPost [MHImagePar]", "MImagePar", "FillImgParPost");
944 MFillH fill7b("MHNewParOnPost [MHNewImagePar]", "MNewImagePar", "FillNewParPost");
945 //MFillH fill9b("MHEffOnTime [MHEffectiveOnTime]", "MTime", "FillEffOnTime");
946 fill1b.SetNameTab("PreCut");
947 fill2b.SetNameTab("PostCut");
948 fill3b.SetNameTab("VsSize");
949 fill3d.SetNameTab("CutT");
950 fill4b.SetNameTab("HilExt");
951 fill5b.SetNameTab("HilSrc");
952 fill6b.SetNameTab("ImgPar");
953 fill7b.SetNameTab("NewPar");
954 //fill9b.SetNameTab("EffOnT");
955 fill1b.SetDrawOption(set.HasOffSequences()||set.IsWobbleMode()?"same":"");
956 fill2b.SetDrawOption(set.HasOffSequences()||set.IsWobbleMode()?"same":"");
957 fill3b.SetDrawOption(set.HasOffSequences()||set.IsWobbleMode()?"same":"");
958 fill3d.SetDrawOption(set.HasOffSequences()||set.IsWobbleMode()?"same":"");
959 fill4b.SetDrawOption(set.HasOffSequences()||set.IsWobbleMode()?"same":"");
960 fill5b.SetDrawOption(set.HasOffSequences()||set.IsWobbleMode()?"same":"");
961 fill6b.SetDrawOption(set.HasOffSequences()||set.IsWobbleMode()?"same":"");
962 fill7b.SetDrawOption(set.HasOffSequences()||set.IsWobbleMode()?"same":"");
963 //fill9b.SetFilter(&fbin);
964
965 fill3d.SetFilter(&ftheta);
966
967 /*
968 MHVsTime hvs("MEffectiveOnTime.fVal");
969 hvs.SetTitle("Effective On-Time vs. Time;;T_{on}");
970 MFillH fillvs(&hvs, "MTimeEffectiveOnTime", "FillOnTime");
971 fillvs.SetNameTab("OnTime");
972 */
973
974 /*
975 MParameterD weight;
976 weight.SetVal(-1);
977 fill2a.SetWeight(&weight);
978 fill3a.SetWeight(&weight);
979 fill4a.SetWeight(&weight);
980 fill5a.SetWeight(&weight);
981 fill6a.SetWeight(&weight);
982 fill7a.SetWeight(&weight);
983 if (fSubstraction)
984 {
985 fill2a.SetNameTab("PostCut-");
986 fill3a.SetNameTab("VsSize-");
987 fill4a.SetNameTab("HilExt-");
988 fill5a.SetNameTab("HilSrc-");
989 fill6a.SetNameTab("ImgPar-");
990 fill7a.SetNameTab("NewPar-");
991 }
992 */
993 MHAlpha *halphaon=CreateNewHist(plist);
994 MFillH falpha2(halphaon, "", "FillHist");
995 MH *hfs=CreateNewHistFS(plist);
996 MFillH ffs2(hfs, "MHillas", "FillFS");
997 MFillH fillphi("MHPhi", "", "FillPhi");
998 fillphi.SetDrawOption("anticut");
999
1000 tlist2.SetNumPasses();
1001
1002 tlist.Replace(&readon);
1003 if (fRndmSrcPos && !set.IsWobbleMode())
1004 tlist2.RemoveFromList(&srcrndm);
1005
1006 MFillH fillsrc(&hsrcpos, "MSrcPosCam", "FillSrcPosCam");
1007 fillsrc.SetNameTab("SrcPos");
1008
1009 if (set.IsWobbleMode()/* && !set.IsMonteCarlo()*/)
1010 tlist2.AddToListBefore(&fillsrc, &hcalc);
1011
1012 MH3 hvs2("MPointingPos.fZd");
1013 hvs2.SetName("Theta");
1014 hvs2.SetTitle("Effective On-Time vs. Zenith Angle;\\Theta [\\circ];T_{on} [s]");
1015
1016 MFillH fillvs2(&hvs2, "", "FillOnTime");
1017 if (!set.IsMonteCarlo())
1018 fillvs2.SetWeight("MEffectiveOnTime");
1019 fillvs2.SetNameTab("OnTime");
1020 fillvs2.SetDrawOption(set.HasOffSequences()||set.IsWobbleMode()?"same":"");
1021
1022 if (!fWriteOnly)
1023 {
1024 tlist.Replace(&fill1b);
1025
1026 tlist2.Replace(&fill2b);
1027 if (fFullDisplay)
1028 {
1029 tlist2.Replace(&fill3b);
1030 tlist2.Replace(&fill3d);
1031 tlist2.Replace(&fill4b);
1032 tlist2.Replace(&fill5b);
1033 tlist2.Replace(&fill6b);
1034 tlist2.Replace(&fill7b);
1035 }
1036 tlist2.Replace(&falpha2);
1037 //if (!fIsMonteCarlo)
1038 // tlist2.Replace(&fill9b);
1039 if (!set.IsWobbleMode() || !fNameHist.IsNull())
1040 tlist2.Replace(&ffs2);
1041 if (set.IsWobbleMode())
1042 {
1043 tlist2.AddToListAfter(&fillphi, &falpha2);
1044 if (!fNameHist.IsNull())
1045 tlist2.RemoveFromList(&ffs);
1046 }
1047
1048 if (!set.IsMonteCarlo())
1049 tlist.Replace(&fillvs2);
1050 else
1051 tlist2.Replace(&fillvs2);
1052 }
1053
1054 // by setting it here it is distributed to all consecutive tasks
1055 tlist.SetAccelerator(MTask::kAccDontReset|MTask::kAccDontTime);
1056
1057 par.SetVal(1);
1058
1059 // Execute first analysis
1060 if (!evtloop.Eventloop(fMaxEvents))
1061 {
1062 *fLog << err << GetDescriptor() << ": Processing of on-sequences failed." << endl;
1063 return -4;
1064 }
1065
1066 // FIXME: Perform fit and plot energy dependant alpha plots
1067 // and fit result to new tabs!
1068 if (!WriteResult(plist, set))
1069 return kFALSE;
1070
1071 *fLog << all << GetDescriptor() << ": Done." << endl;
1072 *fLog << endl << endl;
1073
1074 return kTRUE;
1075}
Note: See TracBrowser for help on using the repository browser.