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

Last change on this file since 6994 was 6994, checked in by tbretz, 19 years ago
*** empty log message ***
File size: 14.0 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.", kTRUE);
194 readreal.AddTree("Drive");
195 //read.AddTree("Trigger");
196 //read.AddTree("Camera");
197 //read.AddTree("CC");
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
263 //MContinue cont1("MEventRate2.fRate/MEventRate.fRate>1.1");
264 //MContinue cont2("MEventRate.fRate/MEventRate2.fRate>1.1");
265
266 MFillH fillvs(&hvs, "MTime", "FillEventRate10s");
267
268 MFillH fill0a(&evt0a, "MSignalCam", "FillSignalCam");
269 MFillH fill0b(&evt0b, "MSignalCam", "FillCntUsedPixels");
270 //MFillH fill0r(&evt0r, "MCerPhotEvt", "FillCntUsedRotated");
271 MFillH fill1("MHHillas", "MHillas", "FillHillas");
272 MFillH fill2("MHHillasExt", "", "FillHillasExt");
273 MFillH fill3("MHHillasSrc", "MHillasSrc", "FillHillasSrc");
274 MFillH fill4("MHImagePar", "MImagePar", "FillImagePar");
275 MFillH fill5("MHNewImagePar", "MNewImagePar", "FillNewImagePar");
276 //MFillH fill6("MHImageParTime","MImageParTime","FillImageParTime");
277 //MFillH fill7("MHNewImagePar2","MNewImagePar2","FillNewImagePar2");
278 MFillH fill8(&h1, "", "FillEventRate");
279 MFillH fill9("MHEffectiveOnTime", "MTime", "FillEffOnTime");
280 //MFillH fillb(&h12, "", "FillEvtRate2");
281 //MFillH fill9("MHCerPhot");
282
283 //fill0r.SetDrawOption("colz");
284 fill8.SetNameTab("EvtRate");
285 fill9.SetNameTab("EffOnTime");
286
287 // ------------------ Setup write task ----------------
288
289 const TString rule(Form("%s{s/_Y_/_I_}", fPathOut.Data()));
290 MWriteRootFile write( 2, rule, fOverwrite?"RECREATE":"NEW");
291 MWriteRootFile writem(2, rule, fOverwrite?"RECREATE":"NEW");
292 // Data
293 write.AddContainer("MHillas", "Events");
294 write.AddContainer("MHillasExt", "Events");
295 write.AddContainer("MHillasSrc", "Events");
296 write.AddContainer("MImagePar", "Events");
297 write.AddContainer("MNewImagePar", "Events");
298 //write.AddContainer("MNewImagePar2", "Events");
299 //write.AddContainer("MImageParTime", "Events");
300 write.AddContainer("MRawEvtHeader", "Events");
301 write.AddContainer("MPointingPos", "Events");
302
303 // Run Header
304 write.AddContainer("MRawRunHeader", "RunHeaders");
305 write.AddContainer("MBadPixelsCam", "RunHeaders");
306 write.AddContainer("MGeomCam", "RunHeaders");
307 write.AddContainer("MObservatory", "RunHeaders");
308
309 // Muon Setup
310 write.AddContainer("BinningRadius", "RunHeaders");
311 write.AddContainer("BinningArcWidth", "RunHeaders");
312 write.AddContainer("BinningRingBroadening", "RunHeaders");
313 write.AddContainer("BinningSizeVsArcRadius", "RunHeaders");
314 write.AddContainer("MMuonSetup", "RunHeaders");
315
316 if (ismc)
317 {
318 // Monte Carlo Data
319 write.AddContainer("MMcEvt", "Events");
320 write.AddContainer("MMcTrig", "Events");
321 // Monte Carlo Run Headers
322 write.AddContainer("MMcRunHeader", "RunHeaders");
323 write.AddContainer("MMcTrigHeader", "RunHeaders");
324 write.AddContainer("MMcFadcHeader", "RunHeaders");
325 write.AddContainer("MMcConfigRunHeader", "RunHeaders");
326 write.AddContainer("MMcCorsikaRunHeader", "RunHeaders");
327 write.AddCopySource("OriginalMC");
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 MTaskList tlist2("Events");
358 tlist2.AddToList(&apply);
359 if (!ismc)
360 tlist2.AddToList(&cont);
361 tlist2.AddToList(&contsw);
362 if (!ismc)
363 {
364 tlist2.AddToList(&rate);
365 //tlist2.AddToList(&rate1);
366 tlist2.AddToList(&fillvs);
367 //tlist2.AddToList(&cont1);
368 //tlist2.AddToList(&cont2);
369 tlist2.AddToList(&fill8);
370 tlist2.AddToList(&fill9);
371 }
372 //tlist2.AddToList(&fillb);
373 tlist2.AddToList(&clean);
374 tlist2.AddToList(&fill0a);
375 tlist2.AddToList(&fill0b);
376 //tlist2.AddToList(&fill0r);
377 tlist2.AddToList(&poscalc);
378 tlist2.AddToList(&hcalc);
379 tlist2.AddToList(&fill1);
380 tlist2.AddToList(&fill2);
381 tlist2.AddToList(&fill3);
382 tlist2.AddToList(&fill4);
383 tlist2.AddToList(&fill5);
384
385 // ----------------------- Muon Analysis ----------------------
386 MF fmuon1("MHillas.fSize>150", "MuonPreCut");
387 MF fmuon2("(MMuonSearchPar.fRadius>180) && (MMuonSearchPar.fRadius<400) &&"
388 "(MMuonSearchPar.fDeviation<45)", "MuonSearchCut");
389 MF fmuon3("(MMuonCalibPar.fArcPhi>175) && (MMuonSearchPar.fDeviation<37)",
390 "MuonFinalCut");
391 MFDataMember fmuon4("MMuonCalibPar.fArcPhi", '>', -0.5, "MuonWriteCut");
392 writem.SetFilter(&fmuon4);
393
394 MMuonSearchParCalc muscalc;
395 muscalc.SetFilter(&fmuon1);
396
397 MMuonCalibParCalc mcalc;
398 mcalc.SetFilter(&fmuon2);
399
400 MFillH fillmuon("MHSingleMuon", "", "FillMuon");
401 MFillH fillmpar("MHMuonPar", "", "FillMuonPar");
402 fillmuon.SetFilter(&fmuon2);
403 fillmpar.SetFilter(&fmuon3);
404
405 if (fMuonAnalysis)
406 {
407 tlist2.AddToList(&fmuon1);
408 tlist2.AddToList(&muscalc);
409 tlist2.AddToList(&fmuon2);
410 tlist2.AddToList(&fillmuon);
411 tlist2.AddToList(&mcalc);
412 tlist2.AddToList(&fmuon3);
413 tlist2.AddToList(&fillmpar);
414 tlist2.AddToList(&fmuon4);
415 tlist2.AddToList(&writem);
416 }
417 // ------------------------------------------------------------
418
419 //tlist2.AddToList(&fill6);
420 //tlist2.AddToList(&fill7);
421 //tlist2.AddToList(&fill9);
422
423 MPointingPosCalc pcalc;
424
425 tlist.AddToList(ismc ? (MTask*)&readmc : (MTask*)&readreal);
426 tlist.AddToList(&pcalc, ismc ? "Events" : "Drive");
427 tlist.AddToList(&tlist2, "Events");
428 tlist.AddToList(&write);
429
430 // Create and setup the eventloop
431 MEvtLoop evtloop(fName);
432 evtloop.SetParList(&plist);
433 evtloop.SetDisplay(fDisplay);
434 evtloop.SetLogStream(fLog);
435 if (!SetupEnv(evtloop))
436 return kFALSE;
437
438 // Execute first analysis
439 if (!evtloop.Eventloop(fMaxEvents))
440 {
441 *fLog << err << GetDescriptor() << ": Failed." << endl;
442 return kFALSE;
443 }
444
445 tlist.PrintStatistics();
446
447 if (!WriteResult())
448 return kFALSE;
449
450 *fLog << all << GetDescriptor() << ": Done." << endl;
451 *fLog << endl << endl;
452
453 return kTRUE;
454}
Note: See TracBrowser for help on using the repository browser.