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