Context Navigation

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

Visit:

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

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

Download in other formats: