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

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