source: trunk/MagicSoft/Mars/mjobs/MJStar.cc@ 8785

Last change on this file since 8785 was 8658, checked in by tbretz, 17 years ago
*** empty log message ***
File size: 16.6 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/2004 <mailto:tbretz@astro.uni-wuerzburg.de>
19!
20! Copyright: MAGIC Software Development, 2000-2007
21!
22!
23\* ======================================================================== */
24
25/////////////////////////////////////////////////////////////////////////////
26//
27// MJStar
28//
29// Resource file entries are case sensitive!
30//
31/////////////////////////////////////////////////////////////////////////////
32#include "MJStar.h"
33
34#include <TEnv.h>
35#include <TFile.h>
36
37#include "MLog.h"
38#include "MLogManip.h"
39
40#include "MDirIter.h"
41#include "MParList.h"
42#include "MTaskList.h"
43#include "MEvtLoop.h"
44
45#include "MStatusDisplay.h"
46
47#include "MHSectorVsTime.h"
48#include "MHEffectiveOnTime.h"
49#include "MHCamEvent.h"
50#include "MBinning.h"
51
52#include "MReadReports.h"
53#include "MReadMarsFile.h"
54#include "MFDataPhrase.h"
55#include "MFilterList.h"
56#include "MFDataMember.h"
57#include "MFDeltaT.h"
58#include "MFSoftwareTrigger.h"
59#include "MContinue.h"
60#include "MGeomApply.h"
61#include "MEventRateCalc.h"
62#include "MImgCleanStd.h"
63#include "MSrcPosCalc.h"
64#include "MSrcPosCorrect.h"
65#include "MHillasCalc.h"
66#include "MMuonSearchParCalc.h"
67#include "MMuonCalibParCalc.h"
68#include "MFillH.h"
69#include "MWriteRootFile.h"
70
71#include "MMuonSetup.h"
72#include "MObservatory.h"
73#include "MPointingPosCalc.h"
74
75ClassImp(MJStar);
76
77using namespace std;
78
79// --------------------------------------------------------------------------
80//
81// Default constructor.
82//
83// Sets fRuns to 0, fExtractor to NULL, fDataCheck to kFALSE
84//
85MJStar::MJStar(const char *name, const char *title) : fMuonAnalysis(kTRUE)
86{
87 fName = name ? name : "MJStar";
88 fTitle = title ? title : "Standard analysis and reconstruction";
89}
90
91Bool_t MJStar::CheckEnvLocal()
92{
93 DisableMuonAnalysis(!GetEnv("MuonAnalysis", fMuonAnalysis));
94 return kTRUE;
95}
96
97Bool_t MJStar::WriteResult()
98{
99 if (fPathOut.IsNull())
100 {
101 *fLog << inf << "No output path specified via SetPathOut - no output written." << endl;
102 return kTRUE;
103 }
104
105 TObjArray cont;
106 cont.Add(const_cast<TEnv*>(GetEnv()));
107
108 if (fDisplay)
109 cont.Add(fDisplay);
110
111 const TString oname = Form("star%08d.root", fSequence.GetSequence());
112 return WriteContainer(cont, oname, "RECREATE");
113}
114
115Bool_t MJStar::Process()
116{
117 if (!fSequence.IsValid())
118 {
119 *fLog << err << "ERROR - Sequence invalid!" << endl;
120 return kFALSE;
121 }
122
123 //if (!CheckEnv())
124 // return kFALSE;
125
126 CheckEnv();
127
128 // --------------------------------------------------------------------------------
129
130 *fLog << inf;
131 fLog->Separator(GetDescriptor());
132 *fLog << "Calculate image parameters of sequence ";
133 *fLog << fSequence.GetFileName() << endl;
134 *fLog << endl;
135
136 // --------------------------------------------------------------------------------
137
138 MDirIter iter;
139 if (fSequence.SetupDatRuns(iter, MSequence::kCalibrated)<=0)
140 return kFALSE;
141
142 // Setup Parlist
143 MParList plist;
144 plist.AddToList(this); // take care of fDisplay!
145
146 MObservatory obs;
147 plist.AddToList(&obs);
148
149 MMuonSetup muonsetup;
150 plist.AddToList(&muonsetup);
151
152 // Setup binnings for muon analysis
153 MBinning bins1("BinningRadius");
154 MBinning bins2("BinningArcWidth");
155 MBinning bins3("BinningRingBroadening");
156 MBinning bins4("BinningSizeVsArcRadius");
157 MBinning bins5("BinningMuonWidth");
158 MBinning bins6("BinningArcPhi");
159 plist.AddToList(&bins1);
160 plist.AddToList(&bins2);
161 plist.AddToList(&bins3);
162 plist.AddToList(&bins4);
163 plist.AddToList(&bins5);
164 plist.AddToList(&bins6);
165
166
167 // Setup Tasklist
168 MTaskList tlist;
169 plist.AddToList(&tlist);
170
171 MReadReports readreal;
172 readreal.AddTree("Events", "MTime.", MReadReports::kMaster);
173 readreal.AddTree("Drive", MReadReports::kRequired);
174 readreal.AddTree("Starguider", MReadReports::kRequired);
175 readreal.AddTree("Currents", MReadReports::kRequired);
176 readreal.AddTree("CC");
177 readreal.AddTree("Trigger");
178 readreal.AddFiles(iter);
179
180 MReadMarsFile readmc("Events");
181 readmc.DisableAutoScheme();
182 readmc.AddFiles(iter);
183
184 // ------------------ Setup general tasks ----------------
185
186 MFDeltaT fdeltat;
187 MContinue cont(&fdeltat, "FilterDeltaT", "Filter events with wrong timing");
188 cont.SetInverted();
189
190 MGeomApply apply; // Only necessary to craete geometry
191 MEventRateCalc rate;
192 rate.SetNumEvents(1200);
193
194 MFSoftwareTrigger swtrig;
195 MContinue contsw(&swtrig, "FilterSwTrigger", "Software trigger");
196 contsw.SetInverted();
197
198 MImgCleanStd clean;
199 clean.SetNamePedPhotCam("MPedPhotFromExtractorRndm");
200
201 MSrcPosCalc poscalc;
202 MHillasCalc hcalc;
203 hcalc.Disable(MHillasCalc::kCalcConc);
204
205 // ------------------ Setup histograms and fill tasks ----------------
206 MHCamEvent evt0a(0, "Cleaned", "Average signal after Cleaning;;S [\\gamma]");
207 MHCamEvent evt0c(0, "Sparkless", "Average signal after Cleaning and Spark cuts;;S [\\gamma]");
208 MHCamEvent evt0d(0, "Sparks", "Average signal after Cleaning for Spark cuts;;S [\\gamma]");
209 MHCamEvent evt0b(0, "UsedPix", "Fraction of Events in which Pixels are used;;Fraction");
210 evt0a.SetErrorSpread(kFALSE);
211 evt0b.SetErrorSpread(kFALSE);
212 evt0c.SetErrorSpread(kFALSE);
213 evt0d.SetErrorSpread(kFALSE);
214 evt0b.SetThreshold(0);
215
216 MFillH fillvs("MHRate", "MTime", "FillEventRate");
217 MFillH fillp1("MHPointing", "MTimeDrive", "FillDrive");
218 MFillH fillp2("MHPointing", "MTimeStarguider", "FillStarguider");
219 fillp1.SetBit(MFillH::kDoNotDisplay);
220 fillp1.SetBit(MFillH::kCanSkip);
221 fillp2.SetBit(MFillH::kCanSkip);
222
223 // Needed in parameter list for ReadEnv
224 MHEffectiveOnTime hontime;
225 plist.AddToList(&hontime);
226
227 MFillH fill0a(&evt0a, "MSignalCam", "FillSignalCam");
228 MFillH fill0b(&evt0b, "MSignalCam", "FillCntUsedPixels");
229 MFillH fill0c(&evt0c, "MSignalCam", "FillSignalCamSparkless");
230 MFillH fill0d(&evt0d, "MSignalCam", "FillSignalCamSparks");
231 MFillH fill1("MHHillas", "MHillas", "FillHillas");
232 MFillH fill2("MHHillasExt", "", "FillHillasExt");
233 MFillH fill3("MHHillasSrc", "MHillasSrc", "FillHillasSrc");
234 MFillH fill4("MHImagePar", "MImagePar", "FillImagePar");
235 MFillH fill5("MHNewImagePar", "MNewImagePar", "FillNewImagePar");
236 MFillH fill9("MHEffectiveOnTime", "MTime", "FillEffOnTime");
237
238 //fillvs.SetNameTab("Rate");
239 fill9.SetNameTab("EffOnTime");
240 fill0c.SetNameTab("Sparkless");
241 fill0d.SetNameTab("Sparks");
242
243 // ------------------ Setup write task ----------------
244
245 // Effective on-time need its own not to be skipped by (eg) image cleaning
246 // Muons needs its own to have a unique SetReadyToSave
247 const TString rule(Form("%s{s/_Y_/_I_}", fPathOut.Data()));
248 MWriteRootFile write( 2, rule, fOverwrite?"RECREATE":"NEW");
249 MWriteRootFile writet(2, rule, fOverwrite?"RECREATE":"NEW"); // EffectiveOnTime
250 MWriteRootFile writem(2, rule, fOverwrite?"RECREATE":"NEW"); // Muons
251 writem.SetName("WriteMuons");
252
253 // Data
254 write.AddContainer("MHillas", "Events");
255 write.AddContainer("MHillasExt", "Events");
256 write.AddContainer("MHillasSrc", "Events");
257 write.AddContainer("MImagePar", "Events");
258 write.AddContainer("MNewImagePar", "Events");
259 write.AddContainer("MRawEvtHeader", "Events");
260 write.AddContainer("MPointingPos", "Events");
261
262 // Run Header
263 write.AddContainer("MRawRunHeader", "RunHeaders");
264// write.AddContainer("MBadPixelsCam", "RunHeaders");
265 write.AddContainer("MGeomCam", "RunHeaders");
266 write.AddContainer("MObservatory", "RunHeaders");
267
268 // Muon Setup
269 write.AddContainer("BinningRadius", "RunHeaders");
270 write.AddContainer("BinningArcWidth", "RunHeaders");
271 write.AddContainer("BinningRingBroadening", "RunHeaders");
272 write.AddContainer("BinningSizeVsArcRadius", "RunHeaders");
273 write.AddContainer("MMuonSetup", "RunHeaders");
274
275 if (fSequence.IsMonteCarlo())
276 {
277 // Monte Carlo Data
278 write.AddContainer("MMcEvt", "Events");
279 write.AddContainer("MMcTrig", "Events");
280 // Monte Carlo Run Headers
281 write.AddContainer("MMcRunHeader", "RunHeaders");
282 write.AddContainer("MMcTrigHeader", "RunHeaders");
283 write.AddContainer("MMcFadcHeader", "RunHeaders");
284 write.AddContainer("MMcConfigRunHeader", "RunHeaders");
285 write.AddContainer("MMcCorsikaRunHeader", "RunHeaders");
286 }
287 else
288 {
289 write.AddContainer("MTime", "Events");
290 // Drive
291 write.AddContainer("MReportDrive", "Drive");
292 write.AddContainer("MTimeDrive", "Drive");
293 // Starguider
294 write.AddContainer("MReportStarguider", "Starguider", kFALSE);
295 write.AddContainer("MTimeStarguider", "Starguider", kFALSE);
296 // Effective On Time
297 writet.AddContainer("MEffectiveOnTime", "EffectiveOnTime");
298 writet.AddContainer("MTimeEffectiveOnTime", "EffectiveOnTime");
299 }
300
301 // What to write in muon tree
302 writem.AddContainer("MMuonSearchPar", "Muons");
303 writem.AddContainer("MMuonCalibPar", "Muons");
304 writem.AddContainer("MHillas", "Muons");
305 writem.AddContainer("MHillasExt", "Muons");
306 writem.AddContainer("MHillasSrc", "Muons");
307 writem.AddContainer("MImagePar", "Muons");
308 writem.AddContainer("MNewImagePar", "Muons");
309 writem.AddContainer("MRawEvtHeader", "Muons");
310 writem.AddContainer("MPointingPos", "Muons");
311 if (fSequence.IsMonteCarlo())
312 {
313 // Monte Carlo Data
314 writem.AddContainer("MMcEvt", "Muons");
315 writem.AddContainer("MMcTrig", "Muons");
316 }
317
318 if (fSequence.IsMonteCarlo())
319 if (fMuonAnalysis)
320 writem.AddCopySource("OriginalMC");
321 else
322 write.AddCopySource("OriginalMC");
323
324 MTaskList tlist2("Events");
325 tlist2.AddToList(&apply);
326 if (!fSequence.IsMonteCarlo())
327 tlist2.AddToList(&cont);
328 tlist2.AddToList(&contsw);
329 if (!fSequence.IsMonteCarlo())
330 {
331 // Calibration events don't enter star at all.
332 tlist2.AddToList(&rate);
333 tlist2.AddToList(&fillvs);
334 tlist2.AddToList(&fill9);
335 tlist2.AddToList(&writet);
336 }
337
338 // Spark cut:
339 // This cut is a little bit different from the default cut in the
340 // ganymed.rc, because the cut in ganymed.rc also suppresses a
341 // little background, while the cut here only shows sparks
342 MFDataPhrase fsparks("log10(MNewImagePar.fConcCOG)<-0.45*(log10(MHillas.fSize)-2.5)-0.24", "SparkCut");
343 //fill0b.SetFilter(&fsparks);
344 fill0c.SetFilter(&fsparks);
345
346 // Inverted spark cut (need not to be a member of the task list
347 // because it fsparks is
348 MFilterList fnsparks(&fsparks);
349 fill0d.SetFilter(&fnsparks);
350
351 tlist2.AddToList(&clean);
352 tlist2.AddToList(&poscalc);
353 tlist2.AddToList(&hcalc);
354 tlist2.AddToList(&fsparks);
355 tlist2.AddToList(&fill0a);
356 tlist2.AddToList(&fill0c);
357 tlist2.AddToList(&fill0d);
358 tlist2.AddToList(&fill0b);
359 tlist2.AddToList(&fill1);
360 tlist2.AddToList(&fill2);
361 tlist2.AddToList(&fill3);
362 tlist2.AddToList(&fill4);
363 tlist2.AddToList(&fill5);
364
365 // ----------------------- Muon Analysis ----------------------
366 // Filter to start muon analysis
367 MFDataPhrase fmuon1("MHillas.fSize>150 && MNewImagePar.fConcCOG<0.1", "MuonPreCut");
368 // Filter to calculate further muon parameters
369 MFDataPhrase fmuon2("(MMuonSearchPar.fRadius>180) && (MMuonSearchPar.fRadius<400) &&"
370 "(MMuonSearchPar.fDeviation<45)", "MuonSearchCut");
371 // Filter to fill the MHMuonPar
372 MFDataPhrase fmuon3("(MMuonCalibPar.fArcPhi>190) && (MMuonSearchPar.fDeviation<35) &&"
373 "(MMuonCalibPar.fArcWidth<0.20) && (MMuonCalibPar.fArcWidth>0.04)",
374 "MuonFinalCut");
375 // Filter to write Muons to Muon tree
376 MFDataMember fmuon4("MMuonCalibPar.fArcPhi", '>', -0.5, "MuonWriteCut");
377 writem.SetFilter(&fmuon4);
378
379 MMuonSearchParCalc muscalc;
380 muscalc.SetFilter(&fmuon1);
381
382 MMuonCalibParCalc mcalc;
383 mcalc.SetFilter(&fmuon2);
384
385 MFillH fillmuon("MHSingleMuon", "", "FillMuon");
386 MFillH fillmpar("MHMuonPar", "", "FillMuonPar");
387 fillmuon.SetFilter(&fmuon2);
388 fillmpar.SetFilter(&fmuon3);
389 fillmuon.SetBit(MFillH::kDoNotDisplay);
390
391 if (fMuonAnalysis)
392 {
393 tlist2.AddToList(&fmuon1);
394 tlist2.AddToList(&muscalc);
395 tlist2.AddToList(&fmuon2);
396 tlist2.AddToList(&fillmuon);
397 tlist2.AddToList(&mcalc);
398 tlist2.AddToList(&fmuon3);
399 tlist2.AddToList(&fillmpar);
400 tlist2.AddToList(&fmuon4);
401 tlist2.AddToList(&writem);
402 }
403
404 // ------------------------------------------------------------
405
406 // Initialize histogram
407 MHSectorVsTime histdc, histrms;
408 histdc.SetNameTime("MTimeCurrents");
409 histdc.SetTitle("Mean of all DC Currents;;<I> [nA]");
410 histdc.SetMinimum(0);
411 histdc.SetMaximum(10);
412 histrms.SetNameTime("MTimeCurrents");
413 histrms.SetTitle("Mean pedestal rms of all pixels;;<\\sigma_{p}> [phe]");
414 histrms.SetType(5);
415 histrms.SetMinimum(0);
416 histrms.SetMaximum(10);
417
418 /*
419 // Define area index [0=inner, 1=outer]
420 // TArrayI inner(1); inner[0] = 0; histdc.SetAreaIndex(inner);
421 */
422
423 // Task to fill the histogram
424 MFillH filldc(&histdc, "MCameraDC", "FillDC");
425 MFillH fillrms(&histrms, "MPedPhotFromExtractorRndm", "FillPedRms");
426 //MFillH filltst("MHTest", "MTime", "FillTest");
427 filldc.SetNameTab("Currents");
428 fillrms.SetNameTab("MeanRms");
429 //filltst.SetNameTab("Test");
430
431 MFillH fillw("MHWeather", "MTimeCC", "FillWeather");
432
433 // instantiate camera histogram containers
434 MHCamEvent evtdt(0, "DT", "Discriminator Threshold;;DT [au]");
435
436 // instantiate fill tasks
437 MFillH filldt(&evtdt, "MCameraTH", "FillDT");
438
439 MHSectorVsTime histipr;
440
441 histipr.SetNameTime("MTimeTrigger");
442 histipr.SetTitle("Mean of all IPR;;<IPR> [Hz]");
443 histipr.SetMinimum(0);
444 //histipr.SetUseMedian();
445
446 // Task to fill the histogram
447 MFillH fillipr(&histipr, "MTriggerIPR", "FillIPR");
448 fillipr.SetNameTab("IPR");
449
450 MPointingPosCalc pcalc;
451
452 // ------------------------------------------------------------
453
454 tlist.AddToList(fSequence.IsMonteCarlo() ? (MTask*)&readmc : (MTask*)&readreal);
455 tlist.AddToList(&pcalc, fSequence.IsMonteCarlo() ? "Events" : "Drive");
456 //tlist.AddToList(&filltst, "Events");
457 tlist.AddToList(&tlist2, "Events");
458 if (!fSequence.IsMonteCarlo())
459 {
460 // initiate task list
461 tlist.AddToList(&fillw, "CC");
462 tlist.AddToList(&fillp1, "Drive");
463 tlist.AddToList(&fillp2, "Starguider");
464 tlist.AddToList(&fillrms, "Currents");
465 tlist.AddToList(&filldc, "Currents");
466 tlist.AddToList(&fillipr, "Trigger");
467 tlist.AddToList(&filldt, "CC");
468 }
469 tlist.AddToList(&write);
470
471 // Create and setup the eventloop
472 MEvtLoop evtloop(fName);
473 evtloop.SetParList(&plist);
474 evtloop.SetDisplay(fDisplay);
475 evtloop.SetLogStream(fLog);
476 if (!SetupEnv(evtloop))
477 return kFALSE;
478
479 // Execute first analysis
480 if (!evtloop.Eventloop(fMaxEvents))
481 {
482 *fLog << err << GetDescriptor() << ": Failed." << endl;
483 return kFALSE;
484 }
485
486 if (!WriteResult())
487 return kFALSE;
488
489 *fLog << all << GetDescriptor() << ": Done." << endl;
490 *fLog << endl << endl;
491
492 return kTRUE;
493}
Note: See TracBrowser for help on using the repository browser.