Context Navigation

source: tags/Mars-V2.1.1/mjobs/MJStar.cc

Visit:

Last change on this file was 9028, checked in by tbretz, 16 years ago
* empty log message *
File size: 19.2 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-2008
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	// Core
38	#include "MLog.h"
39	#include "MLogManip.h"
40
41	#include "MDirIter.h"
42	#include "MParList.h"
43	#include "MTaskList.h"
44	#include "MEvtLoop.h"
45
46	#include "MStatusDisplay.h"
47
48	// Histograms
49	#include "MHSectorVsTime.h"
50	#include "MHEffectiveOnTime.h"
51	#include "MHCamEvent.h"
52	#include "MBinning.h"
53
54	///NEEWWWW
55	#include "MH3.h"
56
57	// Tasks
58	#include "MReadReports.h"
59	#include "MReadMarsFile.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	// Filter
73	#include "MFDataPhrase.h"
74	#include "MFTriggerPattern.h"
75	#include "MFilterList.h"
76	#include "MFDataMember.h"
77	#include "MFDeltaT.h"
78	#include "MFSoftwareTrigger.h"
79
80	// Parameter container
81	#include "MMuonSetup.h"
82	#include "MObservatory.h"
83	//#include "MTriggerPattern.h"
84	#include "MPointingPosCalc.h"
85
86	ClassImp(MJStar);
87
88	using namespace std;
89
90	// --------------------------------------------------------------------------
91	//
92	// Default constructor.
93	//
94	// Sets fRuns to 0, fExtractor to NULL, fDataCheck to kFALSE
95	//
96	MJStar::MJStar(const char name, const char title) : fMuonAnalysis(kTRUE)
97	{
98	fName = name ? name : "MJStar";
99	fTitle = title ? title : "Standard analysis and reconstruction";
100	}
101
102	Bool_t MJStar::CheckEnvLocal()
103	{
104	DisableMuonAnalysis(!GetEnv("MuonAnalysis", fMuonAnalysis));
105	return kTRUE;
106	}
107
108	Bool_t MJStar::WriteResult()
109	{
110	if (fPathOut.IsNull())
111	{
112	*fLog << inf << "No output path specified via SetPathOut - no output written." << endl;
113	return kTRUE;
114	}
115
116	TObjArray cont;
117	cont.Add(const_cast<TEnv*>(GetEnv()));
118	cont.Add(const_cast<MSequence*>(&fSequence));
119
120	if (fDisplay)
121	{
122	TString title = "-- Star: ";
123	title += fSequence.GetSequence();
124	title += " --";
125	fDisplay->SetTitle(title, kFALSE);
126
127	cont.Add(fDisplay);
128	}
129
130	const TString oname = Form("star%08d.root", fSequence.GetSequence());
131	return WriteContainer(cont, oname, "RECREATE");
132	}
133
134	Bool_t MJStar::Process()
135	{
136	if (!fSequence.IsValid())
137	{
138	*fLog << err << "ERROR - Sequence invalid!" << endl;
139	return kFALSE;
140	}
141
142	//if (!CheckEnv())
143	// return kFALSE;
144
145	CheckEnv();
146
147	// --------------------------------------------------------------------------------
148
149	*fLog << inf;
150	fLog->Separator(GetDescriptor());
151	*fLog << "Calculate image parameters of sequence ";
152	*fLog << fSequence.GetFileName() << endl;
153	*fLog << endl;
154
155	// --------------------------------------------------------------------------------
156
157	MDirIter iter;
158	if (fSequence.GetRuns(iter, MSequence::kCalibrated)<=0)
159	return kFALSE;
160
161	// Setup Parlist
162	MParList plist;
163	plist.AddToList(this); // take care of fDisplay!
164
165	MObservatory obs;
166	plist.AddToList(&obs);
167
168	MMuonSetup muonsetup;
169	plist.AddToList(&muonsetup);
170
171	// Setup binnings for muon analysis
172	MBinning bins1("BinningRadius");
173	MBinning bins2("BinningArcWidth");
174	MBinning bins3("BinningRingBroadening");
175	MBinning bins4("BinningSizeVsArcRadius");
176	MBinning bins5("BinningMuonWidth");
177	MBinning bins6("BinningArcPhi");
178	plist.AddToList(&bins1);
179	plist.AddToList(&bins2);
180	plist.AddToList(&bins3);
181	plist.AddToList(&bins4);
182	plist.AddToList(&bins5);
183	plist.AddToList(&bins6);
184
185
186	// Setup Tasklist
187	MTaskList tlist;
188	plist.AddToList(&tlist);
189
190	MReadReports readreal;
191	readreal.AddTree("Events", "MTime.", MReadReports::kMaster);
192	readreal.AddTree("Drive", MReadReports::kRequired);
193	readreal.AddTree("Starguider", MReadReports::kRequired);
194	readreal.AddTree("Camera", MReadReports::kRequired);
195	readreal.AddTree("Currents");
196	readreal.AddTree("CC");
197	readreal.AddTree("Rec");
198	readreal.AddTree("Trigger");
199	readreal.AddTree("Pyrometer");
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", "Average signal after Cleaning;;S [phe]");
229	MHCamEvent evt0c(0, "Sparkless", "Average signal after Cleaning and Spark cuts;;S [phe]");
230	MHCamEvent evt0d(0, "Sparks", "Average signal after Cleaning for Spark cuts;;S [phe]");
231	MHCamEvent evt0b(0, "UsedPix", "Fraction of Events in which Pixels are used;;Fraction");
232	evt0a.SetErrorSpread(kFALSE);
233	evt0b.SetErrorSpread(kFALSE);
234	evt0c.SetErrorSpread(kFALSE);
235	evt0d.SetErrorSpread(kFALSE);
236	evt0b.SetThreshold(0);
237
238	MFillH fillvs("MHRate", "MTime", "FillEventRate");
239	MFillH fillp1("MHPointing", "MTimeDrive", "FillDrive");
240	MFillH fillp2("MHPointing", "MTimeStarguider", "FillStarguider");
241	fillp1.SetBit(MFillH::kDoNotDisplay);
242	fillp1.SetBit(MFillH::kCanSkip);
243	fillp2.SetBit(MFillH::kCanSkip);
244
245	// Needed in parameter list for ReadEnv
246	MHEffectiveOnTime hontime;
247	plist.AddToList(&hontime);
248
249	MFillH fill0a(&evt0a, "MSignalCam", "FillSignalCam");
250	MFillH fill0b(&evt0b, "MSignalCam", "FillCntUsedPixels");
251	MFillH fill0c(&evt0c, "MSignalCam", "FillSignalCamSparkless");
252	MFillH fill0d(&evt0d, "MSignalCam", "FillSignalCamSparks");
253	MFillH fill1("MHHillas", "MHillas", "FillHillas");
254	MFillH fill2("MHHillasExt", "", "FillHillasExt");
255	MFillH fill3("MHHillasSrc", "MHillasSrc", "FillHillasSrc");
256	MFillH fill4("MHImagePar", "MImagePar", "FillImagePar");
257	MFillH fill5("MHNewImagePar", "MNewImagePar", "FillNewImagePar");
258	MFillH fill9("MHEffectiveOnTime", "MTime", "FillEffOnTime");
259
260	//fillvs.SetNameTab("Rate");
261	fill9.SetNameTab("EffOnTime");
262	fill0c.SetNameTab("Sparkless");
263	fill0d.SetNameTab("Sparks");
264
265	// FIXME: This should be replaced by a cut working for all triggers!
266	MTriggerPattern par; // An empty trigger pattern for older files
267	plist.AddToList(&par); // without MTriggerPattern stored in the file
268
269	// For the effective on-time calculation we don't want SUM-only events
270	MFTriggerPattern fsum;
271	fsum.SetDefault(kTRUE);
272	fsum.DenySumTrigger();
273	fsum.AllowTriggerLvl1();
274	//fsum.RequireTriggerLvl1();
275	fill9.SetFilter(&fsum);
276
277	// Plot the trigger pattern rates vs. run-number
278	MH3 hrate("MRawRunHeader.GetFileID", "MTriggerPattern.GetUnprescaled");
279	if (fSequence.IsMonteCarlo())
280	hrate.SetWeight("100./MMcRunHeader.fNumSimulatedShowers");
281	else
282	hrate.SetWeight("1./TMath::Max(MRawRunHeader.GetRunLength,1)");
283	hrate.SetName("Rate");
284	hrate.SetTitle("Event rate after cleaning [Hz];File Id;Trigger Type;");
285	hrate.InitLabels(MH3::kLabelsX);
286	hrate.DefaultLabelY("ERROR");
287	hrate.DefineLabelY( 0, "0"); // 0: No pattern
288	hrate.DefineLabelY( 1, "Trig"); // 1: Lvl1
289	hrate.DefineLabelY( 4, "Trig"); // 4: Lvl2
290	hrate.DefineLabelY( 5, "Trig"); // Lvl2+Lvl1
291	hrate.DefineLabelY(32, "Sum"); // 32: Sum
292	hrate.DefineLabelY(33, "Trig"); // Sum+Lvl1
293	hrate.DefineLabelY(36, "Trig"); // Sum+Lvl2
294	hrate.DefineLabelY(37, "Trig"); // Sum+Lvl1+Lvl2
295
296	MFillH frate(&hrate, "", "FillRate");
297	frate.SetDrawOption("box");
298
299	// ------------------ Setup write task ----------------
300
301	// Effective on-time need its own not to be skipped by (eg) image cleaning
302	// Muons needs its own to have a unique SetReadyToSave
303	const TString rule(Form("%s{s/_Y_/_I_}", fPathOut.Data()));
304	MWriteRootFile write( 2, rule, fOverwrite?"RECREATE":"NEW");
305	MWriteRootFile writet(2, rule, fOverwrite?"RECREATE":"NEW"); // EffectiveOnTime
306	MWriteRootFile writem(2, rule, fOverwrite?"RECREATE":"NEW"); // Muons
307	writem.SetName("WriteMuons");
308
309	// Data
310	write.AddContainer("MHillas", "Events");
311	write.AddContainer("MHillasExt", "Events");
312	write.AddContainer("MHillasSrc", "Events");
313	write.AddContainer("MImagePar", "Events");
314	write.AddContainer("MNewImagePar", "Events");
315	write.AddContainer("MRawEvtHeader", "Events");
316	write.AddContainer("MPointingPos", "Events");
317	write.AddContainer("MTriggerPattern", "Events");
318
319	// Run Header
320	write.AddContainer("MRawRunHeader", "RunHeaders");
321	// write.AddContainer("MBadPixelsCam", "RunHeaders");
322	write.AddContainer("MGeomCam", "RunHeaders");
323	write.AddContainer("MObservatory", "RunHeaders");
324
325	// Muon Setup
326	write.AddContainer("BinningRadius", "RunHeaders");
327	write.AddContainer("BinningArcWidth", "RunHeaders");
328	write.AddContainer("BinningRingBroadening", "RunHeaders");
329	write.AddContainer("BinningSizeVsArcRadius", "RunHeaders");
330	write.AddContainer("MMuonSetup", "RunHeaders");
331
332	if (fSequence.IsMonteCarlo())
333	{
334	// Monte Carlo Data
335	write.AddContainer("MMcEvt", "Events");
336	write.AddContainer("MMcTrig", "Events");
337	// Monte Carlo Run Headers
338	write.AddContainer("MMcRunHeader", "RunHeaders");
339	write.AddContainer("MMcTrigHeader", "RunHeaders");
340	write.AddContainer("MMcFadcHeader", "RunHeaders");
341	write.AddContainer("MMcConfigRunHeader", "RunHeaders");
342	write.AddContainer("MMcCorsikaRunHeader", "RunHeaders");
343	}
344	else
345	{
346	// Event time
347	write.AddContainer("MTime", "Events");
348
349	// Slow-control trees with time information
350	write.AddTree("Drive");
351	write.AddTree("Starguider", kFALSE);
352	write.AddTree("Pyrometer", kFALSE);
353	writet.AddContainer("MEffectiveOnTime", "EffectiveOnTime");
354	writet.AddContainer("MTimeEffectiveOnTime", "EffectiveOnTime");
355	}
356
357	// What to write in muon tree
358	writem.AddContainer("MMuonSearchPar", "Muons");
359	writem.AddContainer("MMuonCalibPar", "Muons");
360	writem.AddContainer("MHillas", "Muons");
361	writem.AddContainer("MHillasExt", "Muons");
362	writem.AddContainer("MHillasSrc", "Muons");
363	writem.AddContainer("MImagePar", "Muons");
364	writem.AddContainer("MNewImagePar", "Muons");
365	writem.AddContainer("MRawEvtHeader", "Muons");
366	writem.AddContainer("MPointingPos", "Muons");
367	if (fSequence.IsMonteCarlo())
368	{
369	// Monte Carlo Data
370	writem.AddContainer("MMcEvt", "Muons");
371	writem.AddContainer("MMcTrig", "Muons");
372	}
373
374	if (fSequence.IsMonteCarlo())
375	if (fMuonAnalysis)
376	writem.AddCopySource("OriginalMC");
377	else
378	write.AddCopySource("OriginalMC");
379
380	MTaskList tlist2("Events");
381	tlist2.AddToList(&apply);
382	if (!fSequence.IsMonteCarlo())
383	tlist2.AddToList(&cont);
384	tlist2.AddToList(&contsw);
385	if (!fSequence.IsMonteCarlo())
386	{
387	// Calibration events don't enter star at all.
388	tlist2.AddToList(&rate);
389	tlist2.AddToList(&fillvs);
390	tlist2.AddToList(&fsum);
391	tlist2.AddToList(&fill9);
392	tlist2.AddToList(&writet);
393	}
394
395	// Spark cut:
396	// This cut is a little bit different from the default cut in the
397	// ganymed.rc, because the cut in ganymed.rc also suppresses a
398	// little background, while the cut here only shows sparks
399	MFDataPhrase fsparks("log10(MNewImagePar.fConcCOG)<-0.45*(log10(MHillas.fSize)-2.5)-0.24", "SparkCut");
400	//fill0b.SetFilter(&fsparks);
401	fill0c.SetFilter(&fsparks);
402
403	// Inverted spark cut (need not to be a member of the task list
404	// because it fsparks is
405	MFilterList fnsparks(&fsparks);
406	fill0d.SetFilter(&fnsparks);
407
408	tlist2.AddToList(&clean);
409	tlist2.AddToList(&poscalc);
410
411	MFillH fillsp("MHSrcPosCam", "MSrcPosCam", "FillSrcPosCam");
412	fillsp.SetNameTab("Src");
413	tlist2.AddToList(&fillsp);
414
415
416	tlist2.AddToList(&hcalc);
417	tlist2.AddToList(&fsparks);
418	tlist2.AddToList(&fill0a);
419	tlist2.AddToList(&fill0c);
420	tlist2.AddToList(&fill0d);
421	tlist2.AddToList(&fill0b);
422	tlist2.AddToList(&frate);
423	tlist2.AddToList(&fill1);
424	tlist2.AddToList(&fill2);
425	tlist2.AddToList(&fill3);
426	tlist2.AddToList(&fill4);
427	tlist2.AddToList(&fill5);
428
429	// ----------------------- Muon Analysis ----------------------
430	// Filter to start muon analysis
431	MFDataPhrase fmuon1("MHillas.fSize>150 && MNewImagePar.fConcCOG<0.1", "MuonPreCut");
432	// Filter to calculate further muon parameters
433	MFDataPhrase fmuon2("(MMuonSearchPar.fRadius>180) && (MMuonSearchPar.fRadius<400) &&"
434	"(MMuonSearchPar.fDeviation<45)", "MuonSearchCut");
435	// Filter to fill the MHMuonPar
436	MFDataPhrase fmuon3("(MMuonCalibPar.fArcPhi>190) && (MMuonSearchPar.fDeviation<35) &&"
437	"(MMuonCalibPar.fArcWidth<0.20) && (MMuonCalibPar.fArcWidth>0.04) &&"
438	"MMuonCalibPar.fRelTimeSigma<1.5",
439	"MuonFinalCut");
440	// Filter to write Muons to Muon tree
441	MFDataMember fmuon4("MMuonCalibPar.fArcPhi", '>', -0.5, "MuonWriteCut");
442	writem.SetFilter(&fmuon4);
443
444	MMuonSearchParCalc muscalc;
445	muscalc.SetFilter(&fmuon1);
446
447	MMuonCalibParCalc mcalc;
448	mcalc.SetFilter(&fmuon2);
449
450	MFillH fillmuon("MHSingleMuon", "", "FillMuon");
451	MFillH fillmpar("MHMuonPar", "", "FillMuonPar");
452	fillmuon.SetFilter(&fmuon2);
453	fillmpar.SetFilter(&fmuon3);
454	fillmuon.SetBit(MFillH::kDoNotDisplay);
455
456	if (fMuonAnalysis)
457	{
458	tlist2.AddToList(&fmuon1);
459	tlist2.AddToList(&muscalc);
460	tlist2.AddToList(&fmuon2);
461	tlist2.AddToList(&fillmuon);
462	tlist2.AddToList(&mcalc);
463	tlist2.AddToList(&fmuon3);
464	tlist2.AddToList(&fillmpar);
465	tlist2.AddToList(&fmuon4);
466	tlist2.AddToList(&writem);
467	}
468
469	// ------------------------------------------------------------
470
471	// Initialize histogram
472	MHSectorVsTime histdc, histrms;
473	histdc.SetNameTime("MTimeCamera");
474	histdc.SetTitle("Average DC currents of all pixels vs time;;<I> [nA]");
475	histdc.SetMinimum(0);
476	histdc.SetMaximum(10);
477	histrms.SetNameTime("MTimeCamera");
478	histrms.SetTitle("Average pedestal rms of all pixels vs time;;<\\sigma_{p}> [phe]");
479	histrms.SetType(5);
480	histrms.SetMinimum(0);
481	histrms.SetMaximum(10);
482
483	/*
484	// Define area index [0=inner, 1=outer]
485	// TArrayI inner(1); inner[0] = 0; histdc.SetAreaIndex(inner);
486	*/
487
488	// Task to fill the histogram
489	MFillH filldc(&histdc, "MCameraDC", "FillDC");
490	MFillH fillrms(&histrms, "MPedPhotFromExtractorRndm", "FillPedRms");
491	//MFillH filltst("MHTest", "MTime", "FillTest");
492	filldc.SetNameTab("Currents");
493	fillrms.SetNameTab("MeanRms");
494	//filltst.SetNameTab("Test");
495
496	MFillH fillw1("MHWeather", "MTimeCC", "FillWeather");
497	MFillH fillw2("MHWeather", "MTimePyrometer", "FillPyrometer");
498	fillw2.SetBit(MFillH::kDoNotDisplay);
499	fillw2.SetBit(MFillH::kCanSkip);
500
501	// instantiate camera histogram containers
502	MHCamEvent evtdt(0, "Thresholds", "Average Discriminator Thresholds;;DT [au]");
503	plist.AddToList(&evtdt);
504
505	// instantiate fill tasks
506	MFillH filldt1("Thresholds", "MCameraTH", "FillDT-CC");
507	MFillH filldt2("Thresholds", "MCameraTH", "FillDT-Rec");
508	filldt1.SetNameTab("DT");
509	filldt2.SetBit(MFillH::kDoNotDisplay);
510	filldt2.SetBit(MFillH::kCanSkip);
511
512	MHSectorVsTime histipr;
513
514	histipr.SetNameTime("MTimeTrigger");
515	histipr.SetTitle("Mean of all IPR;;<IPR> [Hz]");
516	histipr.SetMinimum(0);
517	//histipr.SetUseMedian();
518
519	// Task to fill the histogram
520	MFillH fillipr(&histipr, "MTriggerIPR", "FillIPR");
521	fillipr.SetNameTab("IPR");
522
523	MPointingPosCalc pcalc;
524
525	// ------------------------------------------------------------
526
527	tlist.AddToList(fSequence.IsMonteCarlo() ? (MTask)&readmc : (MTask)&readreal);
528	tlist.AddToList(&pcalc, fSequence.IsMonteCarlo() ? "Events" : "Drive");
529	//tlist.AddToList(&filltst, "Events");
530	tlist.AddToList(&tlist2, "Events");
531	if (!fSequence.IsMonteCarlo())
532	{
533	// initiate task list
534	tlist.AddToList(&fillw1, "CC");
535	tlist.AddToList(&fillw2, "Pyrometer");
536	tlist.AddToList(&fillp1, "Drive");
537	tlist.AddToList(&fillp2, "Starguider");
538	tlist.AddToList(&fillrms, "Camera");
539	tlist.AddToList(&filldc, "Camera");
540	tlist.AddToList(&fillipr, "Trigger");
541	tlist.AddToList(&filldt1, "CC"); // Old files: Receiver information in CC-Tree (Mars<=2.0)
542	tlist.AddToList(&filldt2, "Rec"); // New files: Receiver information in Rec-Tree (Mars >2.0)
543	}
544	if (!HasNullOut())
545	tlist.AddToList(&write);
546
547	// Create and setup the eventloop
548	MEvtLoop evtloop(fName);
549	evtloop.SetParList(&plist);
550	evtloop.SetDisplay(fDisplay);
551	evtloop.SetLogStream(fLog);
552	if (!SetupEnv(evtloop))
553	return kFALSE;
554
555	// Execute first analysis
556	if (!evtloop.Eventloop(fMaxEvents))
557	{
558	*fLog << err << GetDescriptor() << ": Failed." << endl;
559	return kFALSE;
560	}
561
562	if (!WriteResult())
563	return kFALSE;
564
565	*fLog << all << GetDescriptor() << ": Done." << endl << endl << endl;;
566
567	return kTRUE;
568	}

Note: See TracBrowser for help on using the repository browser.

Download in other formats: