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

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