source: tags/Mars-V0.9.4.3/mjobs/MJStar.cc

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