Context Navigation

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

Visit:

Last change on this file since 14853 was 14792, checked in by lyard, 12 years ago
Added fits output for MC data
File size: 34.3 KB

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

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

Download in other formats: