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

Last change on this file was 7115, checked in by tbretz, 20 years ago
*** empty log message ***
File size: 14.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-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 "MHCamEvent.h"
50#include "MHCamEventRot.h"
51#include "MBinning.h"
52
53#include "MReadReports.h"
54#include "MReadMarsFile.h"
55#include "MF.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 "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(Bool_t ismc)
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.GetName() << endl;
140 *fLog << endl;
141
142 // --------------------------------------------------------------------------------
143
144 MDirIter iter;
145 const Int_t n0 = fSequence.SetupDatRuns(iter, MSequence::kCalibrated, fPathData);
146 const Int_t n1 = fSequence.GetNumDatRuns();
147 if (n0==0)
148 {
149 *fLog << err << "ERROR - No input files of sequence found!" << endl;
150 return kFALSE;
151 }
152 if (n0!=n1)
153 {
154 *fLog << err << "ERROR - Number of files found (" << n0 << ") doesn't match number of files in sequence (" << n1 << ")" << endl;
155 if (fLog->GetDebugLevel()>4)
156 {
157 *fLog << dbg << "Files which are searched:" << endl;
158 iter.Print();
159 }
160 return kFALSE;
161 }
162
163 // Setup Parlist
164 MParList plist;
165 plist.AddToList(this); // take care of fDisplay!
166
167 MObservatory obs;
168 plist.AddToList(&obs);
169
170 MMuonSetup muonsetup;
171 plist.AddToList(&muonsetup);
172
173 // Setup binnings for muon analysis
174 MBinning bins1("BinningRadius");
175 MBinning bins2("BinningArcWidth");
176 MBinning bins3("BinningRingBroadening");
177 MBinning bins4("BinningSizeVsArcRadius");
178 MBinning bins5("BinningMuonWidth");
179 MBinning bins6("BinningArcPhi");
180 plist.AddToList(&bins1);
181 plist.AddToList(&bins2);
182 plist.AddToList(&bins3);
183 plist.AddToList(&bins4);
184 plist.AddToList(&bins5);
185 plist.AddToList(&bins6);
186
187
188 // Setup Tasklist
189 MTaskList tlist;
190 plist.AddToList(&tlist);
191
192 MReadReports readreal;
193 readreal.AddTree("Events", "MTime.", MReadReports::kMaster);
194 readreal.AddTree("Drive", MReadReports::kRequired);
195 readreal.AddTree("CC");
196 //read.AddTree("Trigger");
197 //read.AddTree("Camera");
198 //read.AddTree("Currents");
199 readreal.AddFiles(iter);
200
201 MReadMarsFile readmc("Events");
202 readmc.DisableAutoScheme();
203 readmc.AddFiles(iter);
204
205 // ------------------ Setup general tasks ----------------
206
207 MFDeltaT fdeltat;
208 MContinue cont(&fdeltat, "FilterDeltaT", "Filter events with wrong timing");
209 cont.SetInverted();
210
211 MGeomApply apply; // Only necessary to craete geometry
212 MEventRateCalc rate;
213 //MEventRateCalc rate1; // 5min
214 rate.SetNumEvents(1200);
215 //rate1.SetNumEvents(60000);
216 //rate1.SetNameEventRate("MEventRate2");
217 //rate1.SetNameTimeRate("MTimeRate2");
218
219 /*
220 MEventRateCalc rate10000;
221 rate10000.SetNameEventRate("MEventRate10000");
222 rate10000.SetNumEvents(10000);
223 */
224 //MBadPixelsMerge merge(&badpix);
225
226 MFSoftwareTrigger swtrig;
227 MContinue contsw(&swtrig, "FilterSwTrigger", "Software trigger");
228 contsw.SetInverted();
229
230
231 MImgCleanStd clean;
232 clean.SetNamePedPhotCam("MPedPhotFromExtractorRndm");
233
234 MSrcPosCalc poscalc;
235 MHillasCalc hcalc;
236 hcalc.Disable(MHillasCalc::kCalcConc);
237
238 // ------------------ Setup histograms and fill tasks ----------------
239 MHCamEvent evt0a(0, "Cleaned", "Signal after Cleaning;;S [\\gamma]");
240 MHCamEvent evt0b(0, "UsedPix", "Fraction of Events in which Pixels are used;;Fraction");
241 evt0b.SetThreshold(0);
242
243 //MHCamEventRot evt0r("UsedRot", "Pixels marked Used (derotated)");
244 //evt0r.SetThreshold(0);
245
246 MH3 h1("MEventRate.fRate");
247 h1.SetName("MHEventRate");
248 h1.SetTitle("Event Rate distribution;R [Hz];Counts");
249 h1.SetLogy();
250/*
251 MH3 h12("MEventRate10000.fRate");
252 h12.SetName("MHEventRate");
253 h12.SetLogy();
254 */
255 MBinning b1("BinningMHEventRate");
256 b1.SetEdges(150, 0, 1500);
257 plist.AddToList(&b1);
258
259 MHVsTime hvs("MEventRate.fRate");
260 hvs.SetTitle("Rate per 500 events;;R [Hz]");
261 hvs.SetNumEvents(500);
262 hvs.SetMinimum(0);
263
264 //MContinue cont1("MEventRate2.fRate/MEventRate.fRate>1.1");
265 //MContinue cont2("MEventRate.fRate/MEventRate2.fRate>1.1");
266
267 MFillH fillvs(&hvs, "MTime", "FillEventRate10s");
268
269 MFillH fill0a(&evt0a, "MSignalCam", "FillSignalCam");
270 MFillH fill0b(&evt0b, "MSignalCam", "FillCntUsedPixels");
271 //MFillH fill0r(&evt0r, "MCerPhotEvt", "FillCntUsedRotated");
272 MFillH fill1("MHHillas", "MHillas", "FillHillas");
273 MFillH fill2("MHHillasExt", "", "FillHillasExt");
274 MFillH fill3("MHHillasSrc", "MHillasSrc", "FillHillasSrc");
275 MFillH fill4("MHImagePar", "MImagePar", "FillImagePar");
276 MFillH fill5("MHNewImagePar", "MNewImagePar", "FillNewImagePar");
277 //MFillH fill6("MHImageParTime","MImageParTime","FillImageParTime");
278 //MFillH fill7("MHNewImagePar2","MNewImagePar2","FillNewImagePar2");
279 MFillH fill8(&h1, "", "FillEventRate");
280 MFillH fill9("MHEffectiveOnTime", "MTime", "FillEffOnTime");
281 //MFillH fillb(&h12, "", "FillEvtRate2");
282 //MFillH fill9("MHCerPhot");
283
284 //fill0r.SetDrawOption("colz");
285 fill8.SetNameTab("EvtRate");
286 fill9.SetNameTab("EffOnTime");
287
288 // ------------------ Setup write task ----------------
289
290 const TString rule(Form("%s{s/_Y_/_I_}", fPathOut.Data()));
291 MWriteRootFile write( 2, rule, fOverwrite?"RECREATE":"NEW");
292 MWriteRootFile writem(2, rule, fOverwrite?"RECREATE":"NEW");
293 // Data
294 write.AddContainer("MHillas", "Events");
295 write.AddContainer("MHillasExt", "Events");
296 write.AddContainer("MHillasSrc", "Events");
297 write.AddContainer("MImagePar", "Events");
298 write.AddContainer("MNewImagePar", "Events");
299 //write.AddContainer("MNewImagePar2", "Events");
300 //write.AddContainer("MImageParTime", "Events");
301 write.AddContainer("MRawEvtHeader", "Events");
302 write.AddContainer("MPointingPos", "Events");
303
304 // Run Header
305 write.AddContainer("MRawRunHeader", "RunHeaders");
306 write.AddContainer("MBadPixelsCam", "RunHeaders");
307 write.AddContainer("MGeomCam", "RunHeaders");
308 write.AddContainer("MObservatory", "RunHeaders");
309
310 // Muon Setup
311 write.AddContainer("BinningRadius", "RunHeaders");
312 write.AddContainer("BinningArcWidth", "RunHeaders");
313 write.AddContainer("BinningRingBroadening", "RunHeaders");
314 write.AddContainer("BinningSizeVsArcRadius", "RunHeaders");
315 write.AddContainer("MMuonSetup", "RunHeaders");
316
317 if (ismc)
318 {
319 // Monte Carlo Data
320 write.AddContainer("MMcEvt", "Events");
321 write.AddContainer("MMcTrig", "Events");
322 // Monte Carlo Run Headers
323 write.AddContainer("MMcRunHeader", "RunHeaders");
324 write.AddContainer("MMcTrigHeader", "RunHeaders");
325 write.AddContainer("MMcFadcHeader", "RunHeaders");
326 write.AddContainer("MMcConfigRunHeader", "RunHeaders");
327 write.AddContainer("MMcCorsikaRunHeader", "RunHeaders");
328 }
329 else
330 {
331 write.AddContainer("MTime", "Events");
332 // Drive
333 write.AddContainer("MReportDrive", "Drive");
334 write.AddContainer("MTimeDrive", "Drive");
335 // Effective On Time
336 write.AddContainer("MEffectiveOnTime", "EffectiveOnTime");
337 write.AddContainer("MTimeEffectiveOnTime", "EffectiveOnTime");
338 }
339
340 // What to write in muon tree
341 writem.AddContainer("MMuonSearchPar", "Muons");
342 writem.AddContainer("MMuonCalibPar", "Muons");
343 writem.AddContainer("MHillas", "Muons");
344 writem.AddContainer("MHillasExt", "Muons");
345 writem.AddContainer("MHillasSrc", "Muons");
346 writem.AddContainer("MImagePar", "Muons");
347 writem.AddContainer("MNewImagePar", "Muons");
348 writem.AddContainer("MRawEvtHeader", "Muons");
349 writem.AddContainer("MPointingPos", "Muons");
350 if (ismc)
351 {
352 // Monte Carlo Data
353 writem.AddContainer("MMcEvt", "Muons");
354 writem.AddContainer("MMcTrig", "Muons");
355 }
356
357 if (ismc)
358 if (fMuonAnalysis)
359 writem.AddCopySource("OriginalMC");
360 else
361 write.AddCopySource("OriginalMC");
362
363 MTaskList tlist2("Events");
364 tlist2.AddToList(&apply);
365 if (!ismc)
366 tlist2.AddToList(&cont);
367 tlist2.AddToList(&contsw);
368 if (!ismc)
369 {
370 // Calibration events don't enter star at all.
371 tlist2.AddToList(&rate);
372 tlist2.AddToList(&fillvs);
373 tlist2.AddToList(&fill8);
374 tlist2.AddToList(&fill9);
375 }
376 //tlist2.AddToList(&fillb);
377 tlist2.AddToList(&clean);
378 tlist2.AddToList(&fill0a);
379 tlist2.AddToList(&fill0b);
380 //tlist2.AddToList(&fill0r);
381 tlist2.AddToList(&poscalc);
382 tlist2.AddToList(&hcalc);
383 tlist2.AddToList(&fill1);
384 tlist2.AddToList(&fill2);
385 tlist2.AddToList(&fill3);
386 tlist2.AddToList(&fill4);
387 tlist2.AddToList(&fill5);
388
389 // ----------------------- Muon Analysis ----------------------
390 MF fmuon1("MHillas.fSize>150", "MuonPreCut");
391 MF fmuon2("(MMuonSearchPar.fRadius>180) && (MMuonSearchPar.fRadius<400) &&"
392 "(MMuonSearchPar.fDeviation<45)", "MuonSearchCut");
393 MF fmuon3("(MMuonCalibPar.fArcPhi>190) && (MMuonSearchPar.fDeviation<35) &&"
394 "(MMuonCalibPar.fArcWidth<0.20) && (MMuonCalibPar.fArcWidth>0.04)",
395 "MuonFinalCut");
396 MFDataMember fmuon4("MMuonCalibPar.fArcPhi", '>', -0.5, "MuonWriteCut");
397 writem.SetFilter(&fmuon4);
398
399 MMuonSearchParCalc muscalc;
400 muscalc.SetFilter(&fmuon1);
401
402 MMuonCalibParCalc mcalc;
403 mcalc.SetFilter(&fmuon2);
404
405 MFillH fillmuon("MHSingleMuon", "", "FillMuon");
406 MFillH fillmpar("MHMuonPar", "", "FillMuonPar");
407 fillmuon.SetFilter(&fmuon2);
408 fillmpar.SetFilter(&fmuon3);
409
410 if (fMuonAnalysis)
411 {
412 tlist2.AddToList(&fmuon1);
413 tlist2.AddToList(&muscalc);
414 tlist2.AddToList(&fmuon2);
415 tlist2.AddToList(&fillmuon);
416 tlist2.AddToList(&mcalc);
417 tlist2.AddToList(&fmuon3);
418 tlist2.AddToList(&fillmpar);
419 tlist2.AddToList(&fmuon4);
420 tlist2.AddToList(&writem);
421 }
422 // ------------------------------------------------------------
423
424 MFillH fillw("MHWeather", "MTimeCC", "FillWeather");
425 //tlist2.AddToList(&fill6);
426 //tlist2.AddToList(&fill7);
427 //tlist2.AddToList(&fill9);
428
429 MPointingPosCalc pcalc;
430
431 tlist.AddToList(ismc ? (MTask*)&readmc : (MTask*)&readreal);
432 tlist.AddToList(&pcalc, ismc ? "Events" : "Drive");
433 tlist.AddToList(&tlist2, "Events");
434 if (!ismc)
435 tlist.AddToList(&fillw, "CC");
436 tlist.AddToList(&write);
437
438 // Create and setup the eventloop
439 MEvtLoop evtloop(fName);
440 evtloop.SetParList(&plist);
441 evtloop.SetDisplay(fDisplay);
442 evtloop.SetLogStream(fLog);
443 if (!SetupEnv(evtloop))
444 return kFALSE;
445
446 // Execute first analysis
447 if (!evtloop.Eventloop(fMaxEvents))
448 {
449 *fLog << err << GetDescriptor() << ": Failed." << endl;
450 return kFALSE;
451 }
452
453 if (!WriteResult())
454 return kFALSE;
455
456 *fLog << all << GetDescriptor() << ": Done." << endl;
457 *fLog << endl << endl;
458
459 return kTRUE;
460}
Note: See TracBrowser for help on using the repository browser.