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

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