Context Navigation

source: branches/MarsGapdTimeJitter/mjobs/MJSimulation.cc@ 18315

Visit:

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

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

Download in other formats: