Context Navigation

source: branches/Mars_MC/mjobs/MJSimulation.cc@ 17048

Visit:

Last change on this file since 17048 was 17048, checked in by ftemme, 11 years ago
Changed the writing of several header keys for the fitsoutput of ceres from hardcoded values to values from a given runheader, the call of SetupHeaderKeys() is commented out, cause the setting of the header keys should be during MWriteFitsFile::ReInit(), cause the Runheader is not initialized at the moment of the previous call
File size: 35.8 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	write3af.AddContainer("MPhotonEvent","Events");
680
681	write3ar.AddContainer("ElectronicNoise", "RunHeaders");
682	write3ar.AddContainer("IntendedPulsePos", "RunHeaders");
683	write3ar.AddContainer("MRawEvtData", "Events");
684	// It doesn't make much sene to write this information
685	// to the file because the units are internal units not
686	// related to the output samples
687	// if (header.IsDataRun() \|\| fForceTrigger)
688	// write3a.AddContainer("TriggerPos", "Events");
689
690	// I: Dedicated file structure
691	write4a.AddContainer("MHillas", "Events");
692	write4a.AddContainer("MHillasSrc", "Events");
693	write4a.AddContainer("MImagePar", "Events");
694	write4a.AddContainer("MNewImagePar", "Events");
695
696	// Basic MC data
697	write1b.AddContainer("MMcEvtBasic", "OriginalMC");
698	write2b.AddContainer("MMcEvtBasic", "OriginalMC");
699	write3b.AddContainer("MMcEvtBasic", "OriginalMC");
700	write4b.AddContainer("MMcEvtBasic", "OriginalMC");
701
702	// -------------------------------------------------------------------
703
704	MGeomApply apply;
705
706	MSimPSF simpsf;
707
708	MSimGeomCam simgeom;
709	simgeom.SetNameGeomCam("GeomCones");
710
711	MSimCalibrationSignal simcal;
712	simcal.SetNameGeomCam("GeomCones");
713
714	// Sky Quality Meter: 21.82M = 2.02e-4cd/m^2
715	// 1cd = 12.566 lm / 4pi sr
716
717	// FIXME: Simulate photons before cones and QE!
718
719	MSimRandomPhotons simnsb;
720	simnsb.SetFreq(5.8, 0.004); // ph/m^2/nm/sr/ns NSB, 4MHz dark count rate
721	simnsb.SetNameGeomCam("GeomCones");
722
723	// FIXME: How do we handle star-light? We need the rate but we also
724	// need to process it through the mirror!
725
726	MSimAPD simapd;
727	simapd.SetNameGeomCam("GeomCones");
728
729	MSimExcessNoise simexcnoise;
730	MSimBundlePhotons simsum;
731	MSimSignalCam simsignal;
732	MSimCamera simcam;
733	MSimTrigger simtrig;
734	MSimReadout simdaq;
735
736	MContinue conttrig("TriggerPos.fVal<0", "ContTrigger");
737
738	MParameterD pulpos("IntendedPulsePos");
739	// FIXME: Set a default which could be 1/3 of the digitization window
740	// pulpos.SetVal(40); // [ns]
741	plist.AddToList(&pulpos);
742
743	MParameterD trigger("TriggerPos");
744	trigger.SetVal(0);
745	plist.AddToList(&trigger);
746
747	// -------------------------------------------------------------------
748
749	// Remove isolated pixels
750	MImgCleanStd clean(0, 0);
751	//clean.SetCleanLvl0(0); // The level above which isolated pixels are kept
752	clean.SetCleanRings(0);
753	clean.SetMethod(MImgCleanStd::kAbsolute);
754
755	//clean.SetNamePedPhotCam("MPedPhotFromExtractorRndm");
756
757	MHillasCalc hcalc;
758	hcalc.Disable(MHillasCalc::kCalcConc);
759
760	MHCamEvent evt0a(/10/3, "Signal", "Average signal (absolute);;S [ph]");
761	MHCamEvent evt0c(/10/3, "MaxSignal", "Maximum signal (absolute);;S [ph]");
762	MHCamEvent evt0d(/11/8, "ArrTm", "Time after first photon;;T [ns]");
763	evt0a.SetErrorSpread(kFALSE);
764	evt0c.SetCollectMax();
765
766	MContinue cut("", "Cut");
767
768	MFillH fillx0a(&evt0a, "MSignalCam", "FillAvgSignal");
769	MFillH fillx0c(&evt0c, "MSignalCam", "FillMaxSignal");
770	MFillH fillx0d(&evt0d, "MSignalCam", "FillArrTm");
771	MFillH fillx1("MHHillas", "MHillas", "FillHillas");
772	MFillH fillx3("MHHillasSrc", "MHillasSrc", "FillHillasSrc");
773	MFillH fillx4("MHNewImagePar", "MNewImagePar", "FillNewImagePar");
774	MFillH fillth("MHThreshold", "", "FillThreshold");
775	MFillH fillca("MHCollectionArea", "", "FillTrigArea");
776
777	fillth.SetNameTab("Threshold");
778	fillca.SetNameTab("TrigArea");
779
780	// -------------------------------------------------------------------
781
782	if (header.IsDataRun())
783	{
784	tasks.AddToList(&read);
785	tasks.AddToList(&precut);
786	tasks.AddToList(&pointing);
787	tasks.AddToList(&simmmcs);
788	if (!fPathOut.IsNull() && !HasNullOut())
789	{
790	//tasks.AddToList(&write1b);
791	tasks.AddToList(&write2b);
792	if (fCamera)
793	tasks.AddToList(&write3b);
794	if (header.IsDataRun())
795	tasks.AddToList(&write4b);
796	}
797	// if (header.IsPointRun())
798	// tasks.AddToList(&stars);
799	if (1)
800	tasks.AddToList(&simatm); // Here because before fillh1
801	tasks.AddToList(&calcangle);
802	tasks.AddToList(&fillh1);
803	tasks.AddToList(&fillG);
804	if (!header.IsPointRun())
805	{
806	tasks.AddToList(&absapd);
807	tasks.AddToList(&absmir);
808	if (0)
809	tasks.AddToList(&simatm); // FASTER?
810	}
811	tasks.AddToList(&reflect);
812	if (!header.IsPointRun())
813	{
814	tasks.AddToList(&fill0);
815	//tasks.AddToList(&fill1);
816	tasks.AddToList(&fill2);
817	tasks.AddToList(&fill3);
818	tasks.AddToList(&fill4);
819	tasks.AddToList(&fillF1);
820	}
821	tasks.AddToList(&cones);
822	tasks.AddToList(&cones2);
823	//if (header.IsPointRun())
824	// tasks.AddToList(&fillP);
825	tasks.AddToList(&cont1);
826	}
827	// -------------------------------
828	tasks.AddToList(&apply);
829	if (header.IsDataRun())
830	{
831	tasks.AddToList(&simpsf);
832	tasks.AddToList(&simgeom);
833	tasks.AddToList(&cont2);
834	if (!header.IsPointRun())
835	{
836	tasks.AddToList(&fillF2);
837	tasks.AddToList(&fillh2);
838	}
839	tasks.AddToList(&cont3);
840	}
841	if (fCamera)
842	{
843	if (header.IsPedestalRun() \|\| header.IsCalibrationRun())
844	tasks.AddToList(&simcal);
845	tasks.AddToList(&simnsb);
846	tasks.AddToList(&simapd);
847	tasks.AddToList(&simexcnoise);
848	}
849	tasks.AddToList(&simsum);
850	if (fCamera)
851	{
852	tasks.AddToList(&simcam);
853	if (header.IsDataRun() \|\| fForceTrigger)
854	tasks.AddToList(&simtrig);
855	tasks.AddToList(&conttrig);
856	tasks.AddToList(&simdaq);
857	}
858	tasks.AddToList(&simsignal); // What do we do if signal<0?
859	if (!fPathOut.IsNull() && !HasNullOut())
860	{
861	//tasks.AddToList(&write1a);
862	if (!header.IsPedestalRun())
863	tasks.AddToList(&write2a);
864	if (fCamera)
865	tasks.AddToList(&write3a);
866	}
867	// -------------------------------
868	if (fCamera)
869	{
870	// FIXME: MHCollectionArea Trigger Area!
871	if (header.IsDataRun())
872	tasks.AddToList(&fillh3);
873	tasks.AddToList(&filltp);
874	}
875	if (header.IsDataRun())
876	{
877	tasks.AddToList(&fillew);
878	tasks.AddToList(&filled);
879	}
880	if (!header.IsPedestalRun())
881	{
882	tasks.AddToList(&fillx0a);
883	tasks.AddToList(&fillx0c);
884	}
885	if (header.IsDataRun())
886	{
887	tasks.AddToList(&clean);
888	tasks.AddToList(&hcalc);
889	tasks.AddToList(&cut);
890	tasks.AddToList(&fillx0d);
891	tasks.AddToList(&fillx1);
892	//tasks.AddToList(&fillx2);
893	tasks.AddToList(&fillx3);
894	tasks.AddToList(&fillx4);
895	if (!HasNullOut())
896	tasks.AddToList(&write4a);
897	//tasks.AddToList(&fillx5);
898
899	tasks.AddToList(&fillh4);
900	tasks.AddToList(&fillth);
901	tasks.AddToList(&fillca);
902	}
903	//-------------------------------------------
904
905	tasks.SetAccelerator(MTask::kAccDontReset\|MTask::kAccDontTime);
906
907	// Create and setup the eventloop
908	MEvtLoop evtloop(fName);
909	evtloop.SetParList(&plist);
910	evtloop.SetDisplay(fDisplay);
911	evtloop.SetLogStream(&gLog);
912	if (!SetupEnv(evtloop))
913	return kFALSE;
914
915	// FIXME: If pedestal file given use the values from this file
916	//-------------------------------------------
917
918	MGeomCam cam = static_cast<MGeomCam>(env2.GetCont());
919
920	if (binstr.IsDefault())
921	binstr.SetEdgesLin(150, -shape.GetWidth(),
922	header.GetFreqSampling()/1000.*header.GetNumSamples()+shape.GetWidth());
923
924	if (binsd.IsDefault() && cam)
925	binsd.SetEdgesLin(100, 0, cam->GetMaxRadius()*cam->GetConvMm2Deg());
926
927	if (binsd0.IsDefault() && cam)
928	binsd0.SetEdgesLin(100, 0, cam->GetMaxRadius()*cam->GetConvMm2Deg());
929
930
931	header.Print();
932
933	// FIXME: Display acceptance curves!
934
935	if (fDisplay)
936	{
937	TCanvas &c = fDisplay->AddTab("Info");
938	c.Divide(2,2);
939
940	c.cd(1);
941	gPad->SetBorderMode(0);
942	gPad->SetFrameBorderMode(0);
943	gPad->SetGridx();
944	gPad->SetGridy();
945	gROOT->SetSelectedPad(0);
946	shape.DrawClone()->SetBit(kCanDelete);
947
948	if (env0.GetCont() && (header.IsDataRun() \|\| header.IsPointRun()))
949	{
950	c.cd(3);
951	gPad->SetBorderMode(0);
952	gPad->SetFrameBorderMode(0);
953	gPad->SetGridx();
954	gPad->SetGridy();
955	gROOT->SetSelectedPad(0);
956	static_cast<MReflector*>(env0.GetCont())->DrawClone("line")->SetBit(kCanDelete);
957	}
958
959	if (fCamera)
960	{
961	if (env1.GetCont())
962	{
963	c.cd(2);
964	gPad->SetBorderMode(0);
965	gPad->SetFrameBorderMode(0);
966	gPad->SetGridx();
967	gPad->SetGridy();
968	gROOT->SetSelectedPad(0);
969	MHCamera c = new MHCamera(static_cast<MGeomCam&>(env1.GetCont()));
970	c->SetStats(kFALSE);
971	c->SetBit(MHCamera::kNoLegend);
972	c->SetBit(kCanDelete);
973	c->Draw();
974	}
975
976	if (cam)
977	{
978	c.cd(4);
979	gPad->SetBorderMode(0);
980	gPad->SetFrameBorderMode(0);
981	gPad->SetGridx();
982	gPad->SetGridy();
983	gROOT->SetSelectedPad(0);
984	MHCamera c = new MHCamera(cam);
985	c->SetStats(kFALSE);
986	c->SetBit(MHCamera::kNoLegend);
987	c->SetBit(kCanDelete);
988	c->Draw();
989	}
990	}
991
992	TCanvas &d = fDisplay->AddTab("Info2");
993	d.Divide(2,3);
994
995	d.cd(1);
996	gPad->SetBorderMode(0);
997	gPad->SetFrameBorderMode(0);
998	gPad->SetGridx();
999	gPad->SetGridy();
1000	gROOT->SetSelectedPad(0);
1001	splinepde.DrawClone()->SetBit(kCanDelete);
1002
1003	d.cd(3);
1004	gPad->SetBorderMode(0);
1005	gPad->SetFrameBorderMode(0);
1006	gPad->SetGridx();
1007	gPad->SetGridy();
1008	gROOT->SetSelectedPad(0);
1009	splinecones2.DrawClone()->SetBit(kCanDelete);
1010
1011	d.cd(5);
1012	gPad->SetBorderMode(0);
1013	gPad->SetFrameBorderMode(0);
1014	gPad->SetGridx();
1015	gPad->SetGridy();
1016	gROOT->SetSelectedPad(0);
1017	splinemirror.DrawClone()->SetBit(kCanDelete);
1018
1019	d.cd(2);
1020	gPad->SetBorderMode(0);
1021	gPad->SetFrameBorderMode(0);
1022	gPad->SetGridx();
1023	gPad->SetGridy();
1024	gROOT->SetSelectedPad(0);
1025	splinecones.DrawClone()->SetBit(kCanDelete);
1026	//splinecones2.DrawClone("same")->SetBit(kCanDelete);
1027
1028	d.cd(6);
1029	gPad->SetBorderMode(0);
1030	gPad->SetFrameBorderMode(0);
1031	gPad->SetGridx();
1032	gPad->SetGridy();
1033	gROOT->SetSelectedPad(0);
1034	MParSpline all = (MParSpline)splinepde.DrawClone();
1035	//all->SetTitle("Combined acceptance");
1036	all->SetBit(kCanDelete);
1037	if (splinemirror.GetSpline())
1038	all->Multiply(*splinemirror.GetSpline());
1039	if (splinecones2.GetSpline())
1040	all->Multiply(*splinecones2.GetSpline());
1041	}
1042
1043	//-------------------------------------------
1044
1045	// Execute first analysis
1046	if (!evtloop.Eventloop(fMaxEvents))
1047	{
1048	gLog << err << GetDescriptor() << ": Failed." << endl;
1049	return kFALSE;
1050	}
1051
1052	//-------------------------------------------
1053	// FIXME: Display Spline in tab
1054	// FIXME: Display Reflector in tab
1055	// FIXME: Display Routing(?) in tab
1056	// FIXME: Display Camera(s) in tab
1057	//-------------------------------------------
1058
1059	if (!WriteResult(plist, seq, header.GetRunNumber()))
1060	return kFALSE;
1061
1062	*fLog << all << GetDescriptor() << ": Done." << endl << endl << endl;;
1063
1064	return kTRUE;
1065	}

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

Download in other formats: