Context Navigation

source: trunk/MagicSoft/Mars/mjobs/MJSimulation.cc@ 9462

Visit:

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

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

Download in other formats: