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

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