Context Navigation

source: trunk/Mars/mjobs/MJSimulation.cc@ 13954

Visit:

Last change on this file since 13954 was 9931, checked in by tbretz, 14 years ago
Some small fixes in MJSimulation
File size: 30.6 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/2009 <mailto:tbretz@astro.uni-wuerzburg.de>
19	!
20	! Copyright: MAGIC Software Development, 2000-2009
21	!
22	!
23	\* ======================================================================== */
24
25	/////////////////////////////////////////////////////////////////////////////
26	//
27	// MJSimulation
28	//
29	//
30	// Force reading a corsika file even if the footer (RUNE-section) is missing
31	// by setting fForceMode to kTRUE or from the resource file by
32	//
33	// ForceMode: Yes
34	//
35	//
36	// To switch off the simulation of the camera electronics, use:
37	//
38	// Camera: Off
39	//
40	// Note, that the border of the camera and the propertied of the cones
41	// are still simulated (simply because it is fast). Furthermore, this
42	// switches off the trigger for the output, i.e. all data which deposits
43	// at least one photon in at least one pixel is written to the output file.
44	//
45	//
46	// In case of a pedestal or calibration run the artificial trigger can
47	// be "switched off" and the cosmics trrigger "switched on" by setting
48	// fForceTrigger to kTRUE or from the resource file by
49	//
50	// ForceTrigger: Yes
51	//
52	//
53	/////////////////////////////////////////////////////////////////////////////
54	#include "MJSimulation.h"
55
56	#include <TEnv.h>
57	#include <TCanvas.h>
58
59	// Core
60	#include "MLog.h"
61	#include "MLogManip.h"
62
63	#include "MArgs.h"
64	#include "MDirIter.h"
65	#include "MParList.h"
66	#include "MTaskList.h"
67	#include "MEvtLoop.h"
68
69	#include "MStatusDisplay.h"
70
71	// Tasks
72	#include "MCorsikaRead.h"
73	#include "MContinue.h"
74	#include "MFillH.h"
75	#include "MGeomApply.h"
76	#include "MParameterCalc.h"
77	#include "MHillasCalc.h"
78	#include "MImgCleanStd.h"
79	#include "MWriteRootFile.h"
80
81	#include "MSimMMCS.h"
82	#include "MSimAbsorption.h"
83	#include "MSimAtmosphere.h"
84	#include "MSimReflector.h"
85	#include "MSimPointingPos.h"
86	#include "MSimPSF.h"
87	#include "MSimGeomCam.h"
88	#include "MSimSignalCam.h"
89	#include "MSimAPD.h"
90	#include "MSimExcessNoise.h"
91	#include "MSimCamera.h"
92	#include "MSimTrigger.h"
93	#include "MSimReadout.h"
94	#include "MSimRandomPhotons.h"
95	#include "MSimBundlePhotons.h"
96	#include "MSimCalibrationSignal.h"
97
98	// Histograms
99	#include "MBinning.h"
100
101	#include "MHn.h"
102	#include "MHCamera.h"
103	#include "MHCamEvent.h"
104	#include "MHPhotonEvent.h"
105
106	// Container
107	#include "MRawRunHeader.h"
108	#include "MParameters.h"
109	#include "MReflector.h"
110	#include "MParEnv.h"
111	#include "MSpline3.h"
112	#include "MParSpline.h"
113	#include "MGeomCam.h"
114
115	#include "MPedestalCam.h"
116	#include "MPedestalPix.h"
117
118	ClassImp(MJSimulation);
119
120	using namespace std;
121
122	// --------------------------------------------------------------------------
123	//
124	// Default constructor.
125	//
126	// Sets fRuns to 0, fExtractor to NULL, fDataCheck to kFALSE
127	//
128	MJSimulation::MJSimulation(const char name, const char title)
129	: fForceMode(kFALSE), fCamera(kTRUE), fForceTrigger(kFALSE),
130	fOperationMode(kModeData), fRunNumber(-1)
131	{
132	fName = name ? name : "MJSimulation";
133	fTitle = title ? title : "Standard analysis and reconstruction";
134	}
135
136	Bool_t MJSimulation::CheckEnvLocal()
137	{
138	fForceMode = GetEnv("ForceMode", fForceMode);
139	fForceTrigger = GetEnv("ForceTrigger", fForceTrigger);
140	fCamera = GetEnv("Camera", fCamera);
141
142	return kTRUE;
143	}
144	/*
145	TString MJSimulation::GetOutFile(const MSequence &seq) const
146	{
147	return seq.IsValid() ? Form("ceres%08d.root", seq.GetSequence()) : "ceres.root";
148	}
149	*/
150
151	Bool_t MJSimulation::WriteResult(const MParList &plist, const MSequence &seq, Int_t num)
152	{
153	if (fPathOut.IsNull())
154	{
155	*fLog << inf << "No output path specified via SetPathOut - no output written." << endl;
156	return kTRUE;
157	}
158
159	TObjArray cont;
160	cont.Add(const_cast<TEnv*>(GetEnv()));
161	if (seq.IsValid())
162	cont.Add(const_cast<MSequence*>(&seq));
163
164	cont.Add(plist.FindObject("PulseShape"));
165	cont.Add(plist.FindObject("Reflector"));
166	cont.Add(plist.FindObject("MGeomCam"));
167	cont.Add(plist.FindObject("GeomCones"));
168
169	TNamed cmdline("CommandLine", fCommandLine.Data());
170	cont.Add(&cmdline);
171
172	if (fDisplay)
173	{
174	TString title = "-- Ceres";
175	if (seq.IsValid())
176	{
177	title += ": ";
178	title += seq.GetSequence();
179	}
180	title += " --";
181	fDisplay->SetTitle("Ceres", kFALSE);
182
183	cont.Add(fDisplay);
184	}
185
186	const TString name = Form("ceres%08d.root", num);
187	return WriteContainer(cont, name, "RECREATE");
188	}
189
190	void MJSimulation::SetupHist(MHn &hist) const
191	{
192	hist.AddHist("MCorsikaEvtHeader.fTotalEnergy");
193	hist.InitName("Energy");
194	hist.InitTitle("Energy;E [GeV]");
195	hist.SetLog(kTRUE, kTRUE, kFALSE);
196
197	hist.AddHist("MPhotonEvent.GetNumExternal");
198	hist.InitName("Size");
199	hist.InitTitle("Size;S [#]");
200	hist.SetLog(kTRUE, kTRUE, kFALSE);
201
202	/*
203	hist.AddHist("-MCorsikaEvtHeader.fX/100","-MCorsikaEvtHeader.fY/100");
204	hist.SetDrawOption("colz");
205	hist.InitName("Impact;Impact;Impact");
206	hist.InitTitle("Impact;West <--> East [m];South <--> North [m]");
207	hist.SetAutoRange();
208	*/
209
210	hist.AddHist("MCorsikaEvtHeader.GetImpact/100");
211	hist.InitName("Impact");
212	hist.InitTitle("Impact;Impact [m]");
213	hist.SetAutoRange();
214
215	hist.AddHist("MCorsikaEvtHeader.fFirstInteractionHeight/100000");
216	hist.InitName("Height");
217	hist.InitTitle("FirstInteractionHeight;h [km]");
218
219	hist.AddHist("(MCorsikaEvtHeader.fAz+MCorsikaRunHeader.fMagneticFieldAz)TMath::RadToDeg()", "MCorsikaEvtHeader.fZdTMath::RadToDeg()");
220	hist.InitName("SkyOrigin;Az;Zd");
221	hist.InitTitle("Sky Origin;Az [\\circ];Zd [\\circ]");
222	hist.SetDrawOption("colz");
223	hist.SetAutoRange();
224
225	hist.AddHist("IncidentAngle.fVal");
226	hist.InitName("ViewCone");
227	hist.InitTitle("Incident Angle;\\alpha [\\circ]");
228	}
229
230	void MJSimulation::SetupCommonFileStructure(MWriteRootFile &write) const
231	{
232	// Common run headers
233	write.AddContainer("MMcCorsikaRunHeader", "RunHeaders", kFALSE);
234	write.AddContainer("MCorsikaRunHeader", "RunHeaders", kFALSE);
235	write.AddContainer("MRawRunHeader", "RunHeaders");
236	write.AddContainer("MGeomCam", "RunHeaders");
237	write.AddContainer("MMcRunHeader", "RunHeaders");
238
239	// Common events
240	write.AddContainer("MCorsikaEvtHeader", "Events", kFALSE);
241	write.AddContainer("MRawEvtHeader", "Events");
242	write.AddContainer("MMcEvt", "Events");
243	write.AddContainer("IncidentAngle", "Events", kFALSE);
244	}
245
246	Bool_t MJSimulation::Process(const MArgs &args, const MSequence &seq)
247	{
248	/*
249	if (!fSequence.IsValid())
250	{
251	*fLog << err << "ERROR - Sequence invalid!" << endl;
252	return kFALSE;
253	}
254	*/
255
256	// if (!HasWritePermission(CombinePath(fPathOut, GetOutFile(seq))))
257	// return kFALSE;
258
259	*fLog << inf;
260	fLog->Separator(GetDescriptor());
261
262	if (!CheckEnv())
263	return kFALSE;
264
265	if (seq.IsValid())
266	*fLog << fSequence.GetFileName() << endl;
267	else
268	*fLog << "Input: " << args.GetNumArguments() << "-files" << endl;
269	*fLog << endl;
270
271	MDirIter iter;
272	if (seq.IsValid() && seq.GetRuns(iter, MSequence::kCorsika)<=0)
273	{
274	*fLog << err << "ERROR - Sequence valid but without files." << endl;
275	return kFALSE;
276	}
277
278	// --------------------------------------------------------------------------------
279
280	// Setup Parlist
281	MParList plist;
282	plist.AddToList(this); // take care of fDisplay!
283
284	// setup TaskList
285	MTaskList tasks;
286	plist.AddToList(&tasks);
287
288	// --------------------------------------------------------------------------------
289
290	// FIXME: Allow empty GeomCones!
291	MParEnv env0("Reflector");
292	MParEnv env1("GeomCones");
293	MParEnv env2("MGeomCam");
294	env0.SetClassName("MReflector");
295	env1.SetClassName("MGeomCam");
296	env2.SetClassName("MGeomCam");
297	plist.AddToList(&env0);
298	plist.AddToList(&env1);
299	plist.AddToList(&env2);
300
301	plist.FindCreateObj("MPedPhotCam", "MPedPhotFromExtractorRndm");
302
303	MParSpline shape("PulseShape");
304	plist.AddToList(&shape);
305
306	// * FIXME * FIXME * FIXME *
307	plist.FindCreateObj("MMcRunHeader");
308
309	MRawRunHeader header;
310	header.SetValidMagicNumber();
311	//header.InitFadcType(3);
312
313	switch (fOperationMode)
314	{
315	case kModeData:
316	header.SetRunType(MRawRunHeader::kRTMonteCarlo\|MRawRunHeader::kRTData);
317	header.SetRunInfo(0, fRunNumber<0 ? 3 : fRunNumber);
318	break;
319
320	case kModePed:
321	header.SetRunType(MRawRunHeader::kRTMonteCarlo\|MRawRunHeader::kRTPedestal);
322	header.SetSourceInfo("Pedestal");
323	header.SetRunInfo(0, fRunNumber<0 ? 1 : fRunNumber);
324	break;
325
326	case kModeCal:
327	header.SetRunType(MRawRunHeader::kRTMonteCarlo\|MRawRunHeader::kRTCalibration);
328	header.SetSourceInfo("Calibration");
329	header.SetRunInfo(0, fRunNumber<0 ? 2 : fRunNumber);
330	break;
331
332	case kModePointRun:
333	header.SetRunType(MRawRunHeader::kRTMonteCarlo\|MRawRunHeader::kRTPointRun);
334	header.SetRunInfo(0, fRunNumber<0 ? 0 : fRunNumber);
335	break;
336	}
337
338	// FIXME: Move to MSimPointingPos, MSimCalibrationSignal
339	// Can we use this as input for MSimPointingPos?
340	header.SetObservation("On", "MonteCarlo");
341	plist.AddToList(&header);
342	// ++++++++ FIXME FIXME FIXME +++++++++++++
343
344	/*
345	MPedestalCam pedcam;
346	pedcam.Init(geomcam.GetNumPixels());
347	for (UInt_t i=0; i<geomcam.GetNumPixels(); i++)
348	pedcam[i].Set(128./header.GetScale(), 22.5/header.GetScale());
349	plist.AddToList(&pedcam);
350	*/
351
352	// -------------------------------------------------------------------
353
354	MCorsikaRead read;
355	read.SetForceMode(fForceMode);
356
357	if (!seq.IsValid())
358	{
359	for (int i=0; i<args.GetNumArguments(); i++)
360	read.AddFile(args.GetArgumentStr(i));
361	}
362	else
363	read.AddFiles(iter);
364
365	MContinue precut("", "PreCut");
366	precut.IsInverted();
367	precut.SetAllowEmpty();
368
369	MSimMMCS simmmcs;
370
371	MParSpline splinepde("PhotonDetectionEfficiency");
372	MParSpline splinemirror("MirrorReflectivity");
373	MParSpline splinecones("ConesAngularAcceptance");
374	MParSpline splinecones2("ConesTransmission");
375	plist.AddToList(&splinepde);
376	plist.AddToList(&splinemirror);
377	plist.AddToList(&splinecones);
378	plist.AddToList(&splinecones2);
379
380	const TString sin2 = "sin(MCorsikaEvtHeader.fZd)sin(MCorsikaRunHeader.fZdMinTMath::DegToRad())";
381	const TString cos2 = "cos(MCorsikaEvtHeader.fZd)cos(MCorsikaRunHeader.fZdMinTMath::DegToRad())";
382	const TString cos = "cos(MCorsikaEvtHeader.fAz-MCorsikaRunHeader.fAzMin*TMath::DegToRad())";
383
384	const TString form = "acos("+sin2+""+cos+"+"+cos2+")TMath::RadToDeg()";
385
386	MParameterCalc calcangle(form, "CalcIncidentAngle");
387	calcangle.SetNameParameter("IncidentAngle");
388
389	MSimAtmosphere simatm;
390	MSimAbsorption absapd("SimPhotonDetectionEfficiency");
391	MSimAbsorption absmir("SimMirrorReflectivity");
392	MSimAbsorption cones("SimConesAngularAcceptance");
393	MSimAbsorption cones2("SimConesTransmission");
394	absapd.SetParName("PhotonDetectionEfficiency");
395	absmir.SetParName("MirrorReflectivity");
396	cones.SetParName("ConesAngularAcceptance");
397	cones.SetUseTheta();
398	cones2.SetParName("ConesTransmission");
399
400	MSimPointingPos pointing;
401
402	MSimReflector reflect;
403	reflect.SetNameGeomCam("GeomCones");
404	reflect.SetNameReflector("Reflector");
405	// MSimStarField stars;
406
407	MContinue cont1("MPhotonEvent.GetNumPhotons<1", "ContEmpty1");
408	MContinue cont2("MPhotonEvent.GetNumPhotons<1", "ContEmpty2");
409	MContinue cont3("MPhotonEvent.GetNumPhotons<2", "ContEmpty3");
410	cont1.SetAllowEmpty();
411	cont2.SetAllowEmpty();
412	cont3.SetAllowEmpty();
413
414	// -------------------------------------------------------------------
415
416	MBinning binse( 120, 1, 1000000, "BinningEnergy", "log");
417	MBinning binsth( 60, 0.9, 900000, "BinningThreshold", "log");
418	MBinning binsee( 36, 0.9, 900000, "BinningEnergyEst", "log");
419	MBinning binss( 100, 1, 10000000, "BinningSize", "log");
420	// MBinning binsi( 100, -500, 500, "BinningImpact");
421	MBinning binsi( 32, 0, 800, "BinningImpact");
422	MBinning binsh( 150, 0, 50, "BinningHeight");
423	MBinning binsaz(720, -360, 360, "BinningAz");
424	MBinning binszd( 70, 0, 70, "BinningZd");
425	MBinning binsvc( 35, 0, 7, "BinningViewCone");
426	MBinning binsel(150, 0, 50, "BinningTotLength");
427	MBinning binsew(150, 0, 15, "BinningMedLength");
428	MBinning binsd("BinningDistC");
429	MBinning binsd0("BinningDist");
430	MBinning binstr("BinningTrigPos");
431
432	plist.AddToList(&binsee);
433	plist.AddToList(&binse);
434	plist.AddToList(&binsth);
435	plist.AddToList(&binss);
436	plist.AddToList(&binsi);
437	plist.AddToList(&binsh);
438	plist.AddToList(&binszd);
439	plist.AddToList(&binsaz);
440	plist.AddToList(&binsvc);
441	plist.AddToList(&binsel);
442	plist.AddToList(&binsew);
443	plist.AddToList(&binstr);
444	plist.AddToList(&binsd);
445	plist.AddToList(&binsd0);
446
447	MHn mhn1, mhn2, mhn3, mhn4;
448	SetupHist(mhn1);
449	SetupHist(mhn2);
450	SetupHist(mhn3);
451	SetupHist(mhn4);
452
453	MH3 mhtp("TriggerPos.fVal-IntendedPulsePos.fVal-PulseShape.GetWidth");
454	mhtp.SetName("TrigPos");
455	mhtp.SetTitle("Trigger position w.r.t. the first photon hitting a detector");
456
457	MH3 mhew("MPhotonStatistics.fLength");
458	mhew.SetName("TotLength");
459	mhew.SetTitle("Time between first and last photon hitting a detector;L [ns]");
460
461	MH3 mhed("MPhotonStatistics.fTimeMedDev");
462	mhed.SetName("MedLength");
463	mhed.SetTitle("Median deviation (1\\sigma);L [ns]");
464
465	MFillH fillh1(&mhn1, "", "FillCorsika");
466	MFillH fillh2(&mhn2, "", "FillH2");
467	MFillH fillh3(&mhn3, "", "FillH3");
468	MFillH fillh4(&mhn4, "", "FillH4");
469	MFillH filltp(&mhtp, "", "FillTriggerPos");
470	MFillH fillew(&mhew, "", "FillEvtWidth");
471	MFillH filled(&mhed, "", "FillMedDev");
472	fillh1.SetNameTab("Corsika", "Distribution as simulated by Corsika");
473	fillh2.SetNameTab("Detectable", "Distribution of events hit the detector");
474	fillh3.SetNameTab("Triggered", "Distribution of triggered events");
475	fillh4.SetNameTab("Cleaned", "Distribution after cleaning and cuts");
476	filltp.SetNameTab("TrigPos", "TriggerPosition w.r.t the first photon");
477	fillew.SetNameTab("EvtWidth", "Time between first and last photon hitting a detector");
478	filled.SetNameTab("MedDev", "Median deviation of arrival time of photons hitting a detector");
479
480	MHPhotonEvent planeG(1, "HPhotonEventGround"); // Get from MaxImpact
481	MHPhotonEvent plane0(2, "HMirrorPlane0"); // Get from MReflector
482	//MHPhotonEvent plane1(2, "HMirrorPlane1");
483	MHPhotonEvent plane2(2, "HMirrorPlane2");
484	MHPhotonEvent plane3(2, "HMirrorPlane3");
485	MHPhotonEvent plane4(2, "HMirrorPlane4");
486	MHPhotonEvent planeF1(6, "HPhotonEventCamera"); // Get from MGeomCam
487	MHPhotonEvent planeF2(header.IsPointRun()?4:6, "HPhotonEventCones"); // Get from MGeomCam
488
489	plist.AddToList(&planeG);
490	plist.AddToList(&plane0);
491	plist.AddToList(&plane2);
492	plist.AddToList(&plane3);
493	plist.AddToList(&plane4);
494	plist.AddToList(&planeF1);
495	plist.AddToList(&planeF2);;
496
497	//MHPSF psf;
498
499	MFillH fillG(&planeG, "MPhotonEvent", "FillGround");
500	MFillH fill0(&plane0, "MirrorPlane0", "FillReflector");
501	//MFillH fill1(&plane1, "MirrorPlane1", "FillCamShadow");
502	MFillH fill2(&plane2, "MirrorPlane2", "FillCandidates");
503	MFillH fill3(&plane3, "MirrorPlane3", "FillReflected");
504	MFillH fill4(&plane4, "MirrorPlane4", "FillFocal");
505	MFillH fillF1(&planeF1, "MPhotonEvent", "FillCamera");
506	MFillH fillF2(&planeF2, "MPhotonEvent", "FillCones");
507	//MFillH fillP(&psf, "MPhotonEvent", "FillPSF");
508
509	fillG.SetNameTab("Ground", "Photon distribution at ground");
510	fill0.SetNameTab("Reflector", "Photon distribution at reflector plane w/o camera shadow");
511	//fill1.SetNameTab("CamShadow", "Photon distribution at reflector plane w/ camera shadow");
512	fill2.SetNameTab("Candidates", "Can hit photon distribution at reflector plane w/ camera shadow");
513	fill3.SetNameTab("Reflected", "Photon distribution after reflector projected to reflector plane");
514	fill4.SetNameTab("Focal", "Photon distribution at focal plane");
515	fillF1.SetNameTab("Camera", "Photon distribution which hit the detector");
516	fillF2.SetNameTab("Cones", "Photon distribution after cones");
517	//fillP.SetNameTab("PSF", "Photon distribution after cones for all mirrors");
518
519	// -------------------------------------------------------------------
520
521	const char *fmt = Form("s/cer([0-9]+)([0-9][0-9][0-9])/%s\\/%08d.$2_%%c_MonteCarlo$1.root/", Esc(fPathOut).Data(), header.GetRunNumber());
522
523	const TString rule1(Form(fmt, 'R'));
524	const TString rule2(Form(fmt, 'Y'));
525	const TString rule4(Form(fmt, 'I'));
526	const TString rule3(Form(fmt, header.GetRunTypeChar()));
527
528	MWriteRootFile write4a( 2, rule4, fOverwrite?"RECREATE":"NEW", "Star file");
529	MWriteRootFile write4b( 2, rule4, fOverwrite?"RECREATE":"NEW", "Star file");
530	MWriteRootFile write3a( 2, rule3, fOverwrite?"RECREATE":"NEW", "Camera file");
531	MWriteRootFile write3b( 2, rule3, fOverwrite?"RECREATE":"NEW", "Camera file");
532	MWriteRootFile write2a( 2, rule2, fOverwrite?"RECREATE":"NEW", "Signal file");
533	MWriteRootFile write2b( 2, rule2, fOverwrite?"RECREATE":"NEW", "Signal file");
534	MWriteRootFile write1a( 2, rule1, fOverwrite?"RECREATE":"NEW", "Reflector file");
535	MWriteRootFile write1b( 2, rule1, fOverwrite?"RECREATE":"NEW", "Reflector file");
536
537	write1a.SetName("WriteRefData");
538	write1b.SetName("WriteRefMC");
539	write2a.SetName("WriteSigData");
540	write2b.SetName("WriteSigMC");
541	write3a.SetName("WriteCamData");
542	write3b.SetName("WriteCamMC");
543	write4a.SetName("WriteImgData");
544	write4b.SetName("WriteImgMC");
545
546	SetupCommonFileStructure(write1a);
547	SetupCommonFileStructure(write2a);
548	SetupCommonFileStructure(write3a);
549	SetupCommonFileStructure(write4a);
550
551	// R: Dedicated file structure
552	write1a.AddContainer("MPhotonEvent", "Events");
553
554	// Y: Dedicated file structure
555	write2a.AddContainer("MPedPhotFromExtractorRndm", "RunHeaders"); // FIXME: Needed for the signal files to be display in MARS
556	write2a.AddContainer("MSignalCam", "Events");
557
558	// D: Dedicated file structure
559	write3a.AddContainer("ElectronicNoise", "RunHeaders");
560	write3a.AddContainer("IntendedPulsePos", "RunHeaders");
561	write3a.AddContainer("MRawEvtData", "Events");
562	// It doesn't make much sene to write this information
563	// to the file because the units are internal units not
564	// related to the output samples
565	// if (header.IsDataRun() \|\| fForceTrigger)
566	// write3a.AddContainer("TriggerPos", "Events");
567
568	// I: Dedicated file structure
569	write4a.AddContainer("MHillas", "Events");
570	write4a.AddContainer("MHillasSrc", "Events");
571	write4a.AddContainer("MImagePar", "Events");
572	write4a.AddContainer("MNewImagePar", "Events");
573
574	// Basic MC data
575	write1b.AddContainer("MMcEvtBasic", "OriginalMC");
576	write2b.AddContainer("MMcEvtBasic", "OriginalMC");
577	write3b.AddContainer("MMcEvtBasic", "OriginalMC");
578	write4b.AddContainer("MMcEvtBasic", "OriginalMC");
579
580	// -------------------------------------------------------------------
581
582	MGeomApply apply;
583
584	MSimPSF simpsf;
585
586	MSimGeomCam simgeom;
587	simgeom.SetNameGeomCam("GeomCones");
588
589	MSimCalibrationSignal simcal;
590	simcal.SetNameGeomCam("GeomCones");
591
592	// Sky Quality Meter: 21.82M = 2.02e-4cd/m^2
593	// 1cd = 12.566 lm / 4pi sr
594
595	// FIXME: Simulate photons before cones and QE!
596
597	MSimRandomPhotons simnsb;
598	simnsb.SetFreq(5.8, 0.004); // ph/m^2/nm/sr/ns NSB, 4MHz dark count rate
599	simnsb.SetNameGeomCam("GeomCones");
600
601	// FIXME: How do we handle star-light? We need the rate but we also
602	// need to process it through the mirror!
603
604	MSimAPD simapd;
605	simapd.SetNameGeomCam("GeomCones");
606
607	MSimExcessNoise simexcnoise;
608	MSimBundlePhotons simsum;
609	MSimSignalCam simsignal;
610	MSimCamera simcam;
611	MSimTrigger simtrig;
612	MSimReadout simdaq;
613
614	MContinue conttrig("TriggerPos.fVal<0", "ContTrigger");
615
616	MParameterD pulpos("IntendedPulsePos");
617	// FIXME: Set a default which could be 1/3 of the digitization window
618	// pulpos.SetVal(40); // [ns]
619	plist.AddToList(&pulpos);
620
621	MParameterD trigger("TriggerPos");
622	trigger.SetVal(0);
623	plist.AddToList(&trigger);
624
625	// -------------------------------------------------------------------
626
627	// Remove isolated pixels
628	MImgCleanStd clean(0, 0);
629	//clean.SetCleanLvl0(0); // The level above which isolated pixels are kept
630	clean.SetCleanRings(0);
631	clean.SetMethod(MImgCleanStd::kAbsolute);
632
633	//clean.SetNamePedPhotCam("MPedPhotFromExtractorRndm");
634
635	MHillasCalc hcalc;
636	hcalc.Disable(MHillasCalc::kCalcConc);
637
638	MHCamEvent evt0a(/10/3, "Signal", "Average signal (absolute);;S [ph]");
639	MHCamEvent evt0c(/10/3, "MaxSignal", "Maximum signal (absolute);;S [ph]");
640	MHCamEvent evt0d(/11/8, "ArrTm", "Time after first photon;;T [ns]");
641	evt0a.SetErrorSpread(kFALSE);
642	evt0c.SetCollectMax();
643
644	MContinue cut("", "Cut");
645
646	MFillH fillx0a(&evt0a, "MSignalCam", "FillAvgSignal");
647	MFillH fillx0c(&evt0c, "MSignalCam", "FillMaxSignal");
648	MFillH fillx0d(&evt0d, "MSignalCam", "FillArrTm");
649	MFillH fillx1("MHHillas", "MHillas", "FillHillas");
650	MFillH fillx3("MHHillasSrc", "MHillasSrc", "FillHillasSrc");
651	MFillH fillx4("MHNewImagePar", "MNewImagePar", "FillNewImagePar");
652	MFillH fillth("MHThreshold", "", "FillThreshold");
653	MFillH fillca("MHCollectionArea", "", "FillTrigArea");
654
655	fillth.SetNameTab("Threshold");
656	fillca.SetNameTab("TrigArea");
657
658	// -------------------------------------------------------------------
659
660	if (header.IsDataRun())
661	{
662	tasks.AddToList(&read);
663	tasks.AddToList(&precut);
664	tasks.AddToList(&pointing);
665	tasks.AddToList(&simmmcs);
666	if (!fPathOut.IsNull() && !HasNullOut())
667	{
668	//tasks.AddToList(&write1b);
669	tasks.AddToList(&write2b);
670	if (fCamera)
671	tasks.AddToList(&write3b);
672	if (header.IsDataRun())
673	tasks.AddToList(&write4b);
674	}
675	// if (header.IsPointRun())
676	// tasks.AddToList(&stars);
677	if (1)
678	tasks.AddToList(&simatm); // Here because before fillh1
679	tasks.AddToList(&calcangle);
680	tasks.AddToList(&fillh1);
681	tasks.AddToList(&fillG);
682	if (!header.IsPointRun())
683	{
684	tasks.AddToList(&absapd);
685	tasks.AddToList(&absmir);
686	if (0)
687	tasks.AddToList(&simatm); // FASTER?
688	}
689	tasks.AddToList(&reflect);
690	if (!header.IsPointRun())
691	{
692	tasks.AddToList(&fill0);
693	//tasks.AddToList(&fill1);
694	tasks.AddToList(&fill2);
695	tasks.AddToList(&fill3);
696	tasks.AddToList(&fill4);
697	tasks.AddToList(&fillF1);
698	}
699	tasks.AddToList(&cones);
700	tasks.AddToList(&cones2);
701	//if (header.IsPointRun())
702	// tasks.AddToList(&fillP);
703	tasks.AddToList(&cont1);
704	}
705	// -------------------------------
706	tasks.AddToList(&apply);
707	if (header.IsDataRun())
708	{
709	tasks.AddToList(&simpsf);
710	tasks.AddToList(&simgeom);
711	tasks.AddToList(&cont2);
712	if (!header.IsPointRun())
713	{
714	tasks.AddToList(&fillF2);
715	tasks.AddToList(&fillh2);
716	}
717	tasks.AddToList(&cont3);
718	}
719	if (fCamera)
720	{
721	if (header.IsPedestalRun() \|\| header.IsCalibrationRun())
722	tasks.AddToList(&simcal);
723	tasks.AddToList(&simnsb);
724	tasks.AddToList(&simapd);
725	tasks.AddToList(&simexcnoise);
726	}
727	tasks.AddToList(&simsum);
728	if (fCamera)
729	{
730	tasks.AddToList(&simcam);
731	if (header.IsDataRun() \|\| fForceTrigger)
732	tasks.AddToList(&simtrig);
733	tasks.AddToList(&conttrig);
734	tasks.AddToList(&simdaq);
735	}
736	tasks.AddToList(&simsignal); // What do we do if signal<0?
737	if (!fPathOut.IsNull() && !HasNullOut())
738	{
739	//tasks.AddToList(&write1a);
740	if (!header.IsPedestalRun())
741	tasks.AddToList(&write2a);
742	if (fCamera)
743	tasks.AddToList(&write3a);
744	}
745	// -------------------------------
746	if (fCamera)
747	{
748	// FIXME: MHCollectionArea Trigger Area!
749	if (header.IsDataRun())
750	tasks.AddToList(&fillh3);
751	tasks.AddToList(&filltp);
752	}
753	if (header.IsDataRun())
754	{
755	tasks.AddToList(&fillew);
756	tasks.AddToList(&filled);
757	}
758	if (!header.IsPedestalRun())
759	{
760	tasks.AddToList(&fillx0a);
761	tasks.AddToList(&fillx0c);
762	}
763	if (header.IsDataRun())
764	{
765	tasks.AddToList(&clean);
766	tasks.AddToList(&hcalc);
767	tasks.AddToList(&cut);
768	tasks.AddToList(&fillx0d);
769	tasks.AddToList(&fillx1);
770	//tasks.AddToList(&fillx2);
771	tasks.AddToList(&fillx3);
772	tasks.AddToList(&fillx4);
773	if (!HasNullOut())
774	tasks.AddToList(&write4a);
775	//tasks.AddToList(&fillx5);
776
777	tasks.AddToList(&fillh4);
778	tasks.AddToList(&fillth);
779	tasks.AddToList(&fillca);
780	}
781	//-------------------------------------------
782
783	tasks.SetAccelerator(MTask::kAccDontReset\|MTask::kAccDontTime);
784
785	// Create and setup the eventloop
786	MEvtLoop evtloop(fName);
787	evtloop.SetParList(&plist);
788	evtloop.SetDisplay(fDisplay);
789	evtloop.SetLogStream(&gLog);
790	if (!SetupEnv(evtloop))
791	return kFALSE;
792
793	// FIXME: If pedestal file given use the values from this file
794	//-------------------------------------------
795
796	MGeomCam cam = static_cast<MGeomCam>(env2.GetCont());
797
798	if (binstr.IsDefault())
799	binstr.SetEdgesLin(150, -shape.GetWidth(),
800	header.GetFreqSampling()/1000.*header.GetNumSamples()+shape.GetWidth());
801
802	if (binsd.IsDefault() && cam)
803	binsd.SetEdgesLin(100, 0, cam->GetMaxRadius()*cam->GetConvMm2Deg());
804
805	if (binsd0.IsDefault() && cam)
806	binsd0.SetEdgesLin(100, 0, cam->GetMaxRadius()*cam->GetConvMm2Deg());
807
808
809	header.Print();
810
811	// FIXME: Display acceptance curves!
812
813	if (fDisplay)
814	{
815	TCanvas &c = fDisplay->AddTab("Info");
816	c.Divide(2,2);
817
818	c.cd(1);
819	gPad->SetBorderMode(0);
820	gPad->SetFrameBorderMode(0);
821	gPad->SetGridx();
822	gPad->SetGridy();
823	gROOT->SetSelectedPad(0);
824	shape.DrawClone()->SetBit(kCanDelete);
825
826	if (env0.GetCont() && (header.IsDataRun() \|\| header.IsPointRun()))
827	{
828	c.cd(3);
829	gPad->SetBorderMode(0);
830	gPad->SetFrameBorderMode(0);
831	gPad->SetGridx();
832	gPad->SetGridy();
833	gROOT->SetSelectedPad(0);
834	static_cast<MReflector*>(env0.GetCont())->DrawClone("line")->SetBit(kCanDelete);
835	}
836
837	if (fCamera)
838	{
839	if (env1.GetCont())
840	{
841	c.cd(2);
842	gPad->SetBorderMode(0);
843	gPad->SetFrameBorderMode(0);
844	gPad->SetGridx();
845	gPad->SetGridy();
846	gROOT->SetSelectedPad(0);
847	MHCamera c = new MHCamera(static_cast<MGeomCam&>(env1.GetCont()));
848	c->SetStats(kFALSE);
849	c->SetBit(MHCamera::kNoLegend);
850	c->SetBit(kCanDelete);
851	c->Draw();
852	}
853
854	if (cam)
855	{
856	c.cd(4);
857	gPad->SetBorderMode(0);
858	gPad->SetFrameBorderMode(0);
859	gPad->SetGridx();
860	gPad->SetGridy();
861	gROOT->SetSelectedPad(0);
862	MHCamera c = new MHCamera(cam);
863	c->SetStats(kFALSE);
864	c->SetBit(MHCamera::kNoLegend);
865	c->SetBit(kCanDelete);
866	c->Draw();
867	}
868	}
869
870	TCanvas &d = fDisplay->AddTab("Info2");
871	d.Divide(2,3);
872
873	d.cd(1);
874	gPad->SetBorderMode(0);
875	gPad->SetFrameBorderMode(0);
876	gPad->SetGridx();
877	gPad->SetGridy();
878	gROOT->SetSelectedPad(0);
879	splinepde.DrawClone()->SetBit(kCanDelete);
880
881	d.cd(3);
882	gPad->SetBorderMode(0);
883	gPad->SetFrameBorderMode(0);
884	gPad->SetGridx();
885	gPad->SetGridy();
886	gROOT->SetSelectedPad(0);
887	splinecones2.DrawClone()->SetBit(kCanDelete);
888
889	d.cd(5);
890	gPad->SetBorderMode(0);
891	gPad->SetFrameBorderMode(0);
892	gPad->SetGridx();
893	gPad->SetGridy();
894	gROOT->SetSelectedPad(0);
895	splinemirror.DrawClone()->SetBit(kCanDelete);
896
897	d.cd(2);
898	gPad->SetBorderMode(0);
899	gPad->SetFrameBorderMode(0);
900	gPad->SetGridx();
901	gPad->SetGridy();
902	gROOT->SetSelectedPad(0);
903	splinecones.DrawClone()->SetBit(kCanDelete);
904	//splinecones2.DrawClone("same")->SetBit(kCanDelete);
905
906	d.cd(6);
907	gPad->SetBorderMode(0);
908	gPad->SetFrameBorderMode(0);
909	gPad->SetGridx();
910	gPad->SetGridy();
911	gROOT->SetSelectedPad(0);
912	MParSpline all = (MParSpline)splinepde.DrawClone();
913	//all->SetTitle("Combined acceptance");
914	all->SetBit(kCanDelete);
915	if (splinemirror.GetSpline())
916	all->Multiply(*splinemirror.GetSpline());
917	if (splinecones2.GetSpline())
918	all->Multiply(*splinecones2.GetSpline());
919	}
920
921	//-------------------------------------------
922
923	// Execute first analysis
924	if (!evtloop.Eventloop(fMaxEvents))
925	{
926	gLog << err << GetDescriptor() << ": Failed." << endl;
927	return kFALSE;
928	}
929
930	//-------------------------------------------
931	// FIXME: Display Spline in tab
932	// FIXME: Display Reflector in tab
933	// FIXME: Display Routing(?) in tab
934	// FIXME: Display Camera(s) in tab
935	//-------------------------------------------
936
937	if (!WriteResult(plist, seq, header.GetRunNumber()))
938	return kFALSE;
939
940	*fLog << all << GetDescriptor() << ": Done." << endl << endl << endl;;
941
942	return kTRUE;
943	}

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

Download in other formats: