Context Navigation

MJSimulation.cc@ 20094

Visit:

Last change on this file since 20094 was 19848, checked in by tbretz, 5 years ago
Make a public access to MCorsikaRead::CalcNumEvents possible. This allows to use it to check whether all files are valid. Added the new option '--check' to ceres for that. Copy '-tel' to the file names. This is more conevenient in reading them.
File size: 49.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	#include <iostream>
59
60	// Core
61	#include "MLog.h"
62	#include "MLogManip.h"
63
64	#include "MArgs.h"
65	#include "MDirIter.h"
66	#include "MParList.h"
67	#include "MTaskList.h"
68	#include "MEvtLoop.h"
69
70	#include "MStatusDisplay.h"
71
72	// Tasks
73	#include "MCorsikaRead.h"
74	#include "MContinue.h"
75	#include "MFillH.h"
76	#include "MGeomApply.h"
77	#include "MParameterCalc.h"
78	#include "MSrcPosCalc.h"
79	#include "MHillasCalc.h"
80	#include "MImgCleanStd.h"
81	#include "MWriteRootFile.h"
82	#include "MWriteFitsFile.h"
83
84	#include "MSimMMCS.h"
85	#include "MSimAbsorption.h"
86	#include "MSimAtmosphere.h"
87	#include "MSimReflector.h"
88	#include "MSimPointingPos.h"
89	#include "MSimPSF.h"
90	#include "MSimGeomCam.h"
91	#include "MSimSignalCam.h"
92	#include "MSimAPD.h"
93	#include "MSimExcessNoise.h"
94	#include "MSimCamera.h"
95	#include "MSimTrigger.h"
96	#include "MSimReadout.h"
97	#include "MSimRandomPhotons.h"
98	#include "MSimBundlePhotons.h"
99	#include "MSimCalibrationSignal.h"
100
101	// Histograms
102	#include "MBinning.h"
103
104	#include "MHn.h"
105	#include "MHCamera.h"
106	#include "MHCamEvent.h"
107	#include "MHPhotonEvent.h"
108
109	// Container
110	#include "MRawRunHeader.h"
111	#include "MParameters.h"
112	#include "MReflector.h"
113	#include "MParEnv.h"
114	#include "MSpline3.h"
115	#include "MParSpline.h"
116	#include "MGeomCam.h"
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	Bool_t MJSimulation::WriteResult(const MParList &plist, const MSequence &seq, Int_t num)
148	{
149	if (fPathOut.IsNull())
150	{
151	*fLog << inf << "No output path specified via SetPathOut - no output written." << endl;
152	return kTRUE;
153	}
154
155	TObjArray cont;
156	cont.Add(const_cast<TEnv*>(GetEnv()));
157	if (seq.IsValid())
158	cont.Add(const_cast<MSequence*>(&seq));
159
160	cont.Add(plist.FindObject("PulseShape"));
161	cont.Add(plist.FindObject("Reflector"));
162	cont.Add(plist.FindObject("MGeomCam"));
163	cont.Add(plist.FindObject("GeomCones"));
164
165	TNamed cmdline("CommandLine", fCommandLine.Data());
166	cont.Add(&cmdline);
167
168	if (fDisplay)
169	{
170	TString title = "-- Ceres";
171	if (seq.IsValid())
172	{
173	title += ": ";
174	title += seq.GetSequence();
175	}
176	title += " --";
177	fDisplay->SetTitle("Ceres", kFALSE);
178
179	cont.Add(fDisplay);
180	}
181
182	const TString name = Form("ceres%08d.root", num);
183	return WriteContainer(cont, name, "RECREATE");
184	}
185
186	void MJSimulation::SetupHist(MHn &hist) const
187	{
188	hist.AddHist("MCorsikaEvtHeader.fTotalEnergy");
189	hist.InitName("Energy");
190	hist.InitTitle("Energy;E [GeV]");
191	hist.SetLog(kTRUE, kTRUE, kFALSE);
192
193	hist.AddHist("MPhotonEvent.GetNumExternal");
194	hist.InitName("Size");
195	hist.InitTitle("Size;S [#]");
196	hist.SetLog(kTRUE, kTRUE, kFALSE);
197
198	/*
199	hist.AddHist("-MCorsikaEvtHeader.fX/100","-MCorsikaEvtHeader.fY/100");
200	hist.SetDrawOption("colz");
201	hist.InitName("Impact;Impact;Impact");
202	hist.InitTitle("Impact;West <--> East [m];South <--> North [m]");
203	hist.SetAutoRange();
204	*/
205
206	hist.AddHist("MCorsikaEvtHeader.GetImpact/100");
207	hist.InitName("Impact");
208	hist.InitTitle("Impact;Impact [m]");
209	hist.SetAutoRange();
210
211	hist.AddHist("MCorsikaEvtHeader.fFirstInteractionHeight/100000");
212	hist.InitName("Height");
213	hist.InitTitle("FirstInteractionHeight;h [km]");
214
215	hist.AddHist("(MCorsikaEvtHeader.fAz+MCorsikaRunHeader.fMagneticFieldAz)TMath::RadToDeg()", "MCorsikaEvtHeader.fZdTMath::RadToDeg()");
216	hist.InitName("SkyOrigin;Az;Zd");
217	hist.InitTitle("Sky Origin;Az [#circ];Zd [#circ]");
218	hist.SetDrawOption("colz");
219	hist.SetAutoRange();
220
221	hist.AddHist("IncidentAngle.fVal");
222	hist.InitName("ViewCone");
223	hist.InitTitle("Incident Angle;#alpha [#circ]");
224	}
225
226	template<class T>
227	void MJSimulation::SetupCommonFileStructure(T &write) const
228	{
229	// Common run headers
230	write.AddContainer("MMcCorsikaRunHeader", "RunHeaders", kFALSE, 1);
231	write.AddContainer("MCorsikaRunHeader", "RunHeaders", kFALSE, 1);
232	write.AddContainer("MRawRunHeader", "RunHeaders", kTRUE, 1);
233	write.AddContainer("MGeomCam", "RunHeaders", kTRUE, 1);
234	write.AddContainer("MMcRunHeader", "RunHeaders", kTRUE, 1);
235
236	// Common events
237	write.AddContainer("MCorsikaEvtHeader", "Events", kFALSE);
238	write.AddContainer("MRawEvtHeader", "Events");
239	write.AddContainer("IncidentAngle", "Events", kFALSE);
240	write.AddContainer("MPointingPos", "Events", kFALSE);
241	write.AddContainer("MSimSourcePos", "Events", kFALSE);
242	write.AddContainer("MSrcPosCam", "Events", kFALSE);
243
244	if (fOperationMode!=kModeData)
245	return;
246
247	write.AddContainer("MMcEvt", "Events", kFALSE);
248	write.AddContainer("MMcEvtBasic", "Events", kFALSE);
249	}
250
251	void MJSimulation::SetupHeaderKeys(MWriteFitsFile &write,MRawRunHeader &header) const
252	{
253	const MTime now(-1);
254	write.SetHeaderKey("ISMC",true,"Bool if File is Montecarlo File");
255	write.SetHeaderKey("TELESCOP", "FACT", "");
256	write.SetHeaderKey("PACKAGE", "MARS Cheobs", "");
257	write.SetHeaderKey("VERSION", "1.0", "");
258	write.SetHeaderKey("CREATOR", "Ceres", "");
259	write.SetHeaderKey("EXTREL", 1., "");
260	write.SetHeaderKey("COMPILED", __DATE__" " __TIME__, "");
261	write.SetHeaderKey("REVISION", "0", "");
262	write.SetHeaderKey("ORIGIN", "FACT", "");
263	write.SetHeaderKey("DATE", now.GetStringFmt("%Y-%m-%dT%H:%M:%S").Data(), "");
264	write.SetHeaderKey("NIGHT", now.GetNightAsInt(), "");
265	write.SetHeaderKey("TIMESYS", "UTC", "");
266	write.SetHeaderKey("TIMEUNIT", "d", "");
267	write.SetHeaderKey("MJDREF", 40587, "");
268	//write.SetHeaderKey("BLDVER", 1, "");
269	write.SetHeaderKey("RUNID", header.GetRunNumber(), "");
270	write.SetHeaderKey("NBOARD", 40, "");
271	write.SetHeaderKey("NPIX", header.GetNumPixel(), "");
272	write.SetHeaderKey("NROI", header.GetNumSamplesHiGain(), "");
273	write.SetHeaderKey("NROITM", 0, "");
274	write.SetHeaderKey("TMSHIFT", 0, "");
275	write.SetHeaderKey("CAMERA", "MGeomCamFACT", "");
276	write.SetHeaderKey("DAQ", "DRS4", "");
277
278	// FTemme: ADCRANGE and ADC have to be calculated, using the values for
279	// the fadctype.
280	// write.SetHeaderKey("ADCRANGE", 2000, "Dynamic range in mV");
281	// write.SetHeaderKey("ADC", 12, "Resolution in bits");
282
283	switch(header.GetRunType())
284	{
285	case MRawRunHeader::kRTData\|MRawRunHeader::kRTMonteCarlo:
286	write.SetHeaderKey("RUNTYPE", "data", "");
287	break;
288	case MRawRunHeader::kRTPedestal\|MRawRunHeader::kRTMonteCarlo:
289	write.SetHeaderKey("RUNTYPE", "pedestal", "");
290	break;
291	case MRawRunHeader::kRTCalibration\|MRawRunHeader::kRTMonteCarlo:
292	write.SetHeaderKey("RUNTYPE", "calibration", "");
293	break;
294	}
295	// write.SetHeaderKey("ID", 777, "Board 0: Board ID");
296	// write.SetHeaderKey("FMVER", 532, "Board 0: Firmware Version");
297	// write.SetHeaderKey("DNA", "0", "");
298	// write.SetHeaderKey("BOARD", 0, "");
299	// write.SetHeaderKey("PRESC", 40, "");
300	// write.SetHeaderKey("PHASE", 0, "");
301	// write.SetHeaderKey("DAC0", 26500, "");
302	// write.SetHeaderKey("DAC1", 0, "");
303	// write.SetHeaderKey("DAC2", 0, "");
304	// write.SetHeaderKey("DAC3", 0, "");
305	// write.SetHeaderKey("DAC4", 28800, "");
306	// write.SetHeaderKey("DAC5", 28800, "");
307	// write.SetHeaderKey("DAC6", 28800, "");
308	// write.SetHeaderKey("DAC7", 28800, "");
309	write.SetHeaderKey("REFCLK", header.GetFreqSampling(), "");
310	write.SetHeaderKey("DRSCALIB", false, "");
311	// write.SetHeaderKey("TSTARTI", 0, "");
312	// write.SetHeaderKey("TSTARTF", 0., "");
313	// write.SetHeaderKey("TSTOPI", 0, "");
314	// write.SetHeaderKey("TSTOPF", 0., "");
315	// write.SetHeaderKey("DATE-OBS", "1970-01-01T00:00:00", "");
316	// write.SetHeaderKey("DATE-END", "1970-01-01T00:00:00", "");
317	// write.SetHeaderKey("NTRG", 0, "");
318	// write.SetHeaderKey("NTRGPED", 0, "");
319	// write.SetHeaderKey("NTRGLPE", 0, "");
320	// write.SetHeaderKey("NTRGTIM", 0, "");
321	// write.SetHeaderKey("NTRGLPI", 0, "");
322	// write.SetHeaderKey("NTRGEXT1", 0, "");
323	// write.SetHeaderKey("NTRGEXT2", 0, "");
324	// write.SetHeaderKey("NTRGMISC", 0, "");
325	}
326
327	void MJSimulation::SetupVetoColumns(MWriteFitsFile &write) const
328	{
329	write.VetoColumn("MParameterD.fVal");
330	write.VetoColumn("MRawEvtData.fLoGainPixId");
331	write.VetoColumn("MRawEvtData.fLoGainFadcSamples");
332	write.VetoColumn("MRawEvtData.fABFlags");
333	write.VetoColumn("MRawEvtHeader.fNumTrigLvl2");
334	//write.VetoColumn("MRawEvtHeader.fTrigPattern");
335	write.VetoColumn("MRawEvtHeader.fNumLoGainOn");
336	}
337
338	// Setup the header accordingly to the used operation mode
339	void MJSimulation::SetupHeaderOperationMode(MRawRunHeader &header) const
340	{
341	switch (fOperationMode)
342	{
343	case kModeData:
344	header.SetRunType(MRawRunHeader::kRTMonteCarlo\|MRawRunHeader::kRTData);
345	header.SetRunInfo(0, fRunNumber<0 ? 3 : fRunNumber);
346	break;
347
348	case kModePed:
349	header.SetRunType(MRawRunHeader::kRTMonteCarlo\|MRawRunHeader::kRTPedestal);
350	header.SetSourceInfo("Pedestal");
351	header.SetRunInfo(0, fRunNumber<0 ? 1 : fRunNumber);
352	break;
353
354	case kModeCal:
355	header.SetRunType(MRawRunHeader::kRTMonteCarlo\|MRawRunHeader::kRTCalibration);
356	header.SetSourceInfo("Calibration");
357	header.SetRunInfo(0, fRunNumber<0 ? 2 : fRunNumber);
358	break;
359
360	case kModePointRun:
361	header.SetRunType(MRawRunHeader::kRTMonteCarlo\|MRawRunHeader::kRTPointRun);
362	header.SetRunInfo(0, fRunNumber<0 ? 0 : fRunNumber);
363	break;
364	}
365	}
366
367	void MJSimulation::CreateBinningObjects(MParList &plist) const
368	{
369	MBinning binse = (MBinning) plist.FindCreateObj("MBinning","BinningEnergy");
370	binse->SetEdges(120, 1, 1000000,"log");
371	MBinning binsth = (MBinning) plist.FindCreateObj("MBinning","BinningThreshold");
372	binsth->SetEdges( 60, 0.9, 900000, "log");
373	MBinning binsee = (MBinning) plist.FindCreateObj("MBinning","BinningEnergyEst");
374	binsee->SetEdges( 36, 0.9, 900000, "log");
375	MBinning binss = (MBinning) plist.FindCreateObj("MBinning","BinningSize");
376	binss->SetEdges( 100, 1, 10000000, "log");
377	MBinning binsi = (MBinning) plist.FindCreateObj("MBinning","BinningImpact");
378	binsi->SetEdges( 32, 0, 800);
379	MBinning binsh = (MBinning) plist.FindCreateObj("MBinning","BinningHeight");
380	binsh->SetEdges( 150, 0, 50);
381	MBinning binsaz = (MBinning) plist.FindCreateObj("MBinning","BinningAz");
382	binsaz->SetEdges(720, -360, 360);
383	MBinning binszd = (MBinning) plist.FindCreateObj("MBinning","BinningZd");
384	binszd->SetEdges( 70, 0, 70);
385	MBinning binsvc = (MBinning) plist.FindCreateObj("MBinning","BinningViewCone");
386	binsvc->SetEdges( 35, 0, 7);
387	MBinning binsel = (MBinning) plist.FindCreateObj("MBinning","BinningTotLength");
388	binsel->SetEdges(150, 0, 50);
389	MBinning binsew = (MBinning) plist.FindCreateObj("MBinning","BinningMedLength");
390	binsew->SetEdges(150, 0, 15);
391	plist.FindCreateObj("MBinning","BinningWidth");
392	plist.FindCreateObj("MBinning","BinningLength");
393	plist.FindCreateObj("MBinning","BinningDistC");
394	plist.FindCreateObj("MBinning","BinningDist");
395	plist.FindCreateObj("MBinning","BinningTrigPos");
396	}
397
398	Bool_t MJSimulation::Process(const MArgs &args, const MSequence &seq, Bool_t checkonly)
399	{
400	// --------------------------------------------------------------------------------
401	// --------------------------------------------------------------------------------
402	// Initialization
403	// --------------------------------------------------------------------------------
404	// --------------------------------------------------------------------------------
405
406	// --------------------------------------------------------------------------------
407	// Logging output:
408	// --------------------------------------------------------------------------------
409	// - description of the job itself
410	// - list of the to processed sequence
411	if (!checkonly)
412	fLog->Separator(GetDescriptor());
413
414	*fLog << inf;
415	if (!CheckEnv())
416	return kFALSE;
417
418	if (seq.IsValid())
419	*fLog << fSequence.GetFileName() << endl;
420	else
421	*fLog << "Input: " << args.GetNumArguments() << " file(s)" << endl;
422	*fLog << endl;
423
424	MDirIter iter;
425	if (seq.IsValid() && seq.GetRuns(iter, MSequence::kCorsika)<=0)
426	{
427	*fLog << err << "ERROR - Sequence valid but without files." << endl;
428	return kFALSE;
429	}
430	// --------------------------------------------------------------------------------
431	// Setup MParList and MTasklist
432	// --------------------------------------------------------------------------------
433	MParList plist;
434	plist.AddToList(this); // take care of fDisplay!
435	// setup TaskList
436	MTaskList tasks;
437	plist.AddToList(&tasks);
438
439	// --------------------------------------------------------------------------------
440	// --------------------------------------------------------------------------------
441	// Parameter Container Setup
442	// --------------------------------------------------------------------------------
443	// --------------------------------------------------------------------------------
444
445	// --------------------------------------------------------------------------------
446	// Setup container for the reflector, the cones and the camera geometry
447	// --------------------------------------------------------------------------------
448	// FIXME: Allow empty GeomCones!
449	MParEnv env0("Reflector");
450	MParEnv env1("GeomCones");
451	MParEnv env2("MGeomCam");
452	env0.SetClassName("MOptics");
453	env1.SetClassName("MGeomCam");
454	env2.SetClassName("MGeomCam");
455	plist.AddToList(&env0);
456	plist.AddToList(&env1);
457	plist.AddToList(&env2);
458	// --------------------------------------------------------------------------------
459	// Setup container for the "camera array" of the extracted number of photons
460	// --------------------------------------------------------------------------------
461	// from ExtractorRndm
462	plist.FindCreateObj("MPedPhotCam", "MPedPhotFromExtractorRndm");
463	// --------------------------------------------------------------------------------
464	// Setup spline containers for:
465	// --------------------------------------------------------------------------------
466	// - the pulse shape
467	// - the different photon acceptance splines
468	MParSpline shape("PulseShape");
469	MParSpline splinepde("PhotonDetectionEfficiency");
470	MParSpline splinemirror("MirrorReflectivity");
471	MParSpline splinecones("ConesAngularAcceptance");
472	MParSpline splinecones2("ConesTransmission");
473	MParSpline splineAdditionalPhotonAcceptance("AdditionalPhotonAcceptance");
474	plist.AddToList(&shape);
475	plist.AddToList(&splinepde);
476	plist.AddToList(&splinemirror);
477	plist.AddToList(&splinecones);
478	plist.AddToList(&splinecones2);
479	plist.AddToList(&splineAdditionalPhotonAcceptance);
480
481	// --------------------------------------------------------------------------------
482	// Setup container for the MC run header and the Raw run header
483	// --------------------------------------------------------------------------------
484	plist.FindCreateObj("MMcRunHeader");
485
486	MRawRunHeader header;
487	header.SetValidMagicNumber();
488	SetupHeaderOperationMode(header);
489	// FIXME: Move to MSimPointingPos, MSimCalibrationSignal
490	// Can we use this as input for MSimPointingPos?
491	header.SetObservation("On", "MonteCarlo");
492	plist.AddToList(&header);
493
494	// --------------------------------------------------------------------------------
495	// Setup container for the intended pulse position and for the trigger position
496	// --------------------------------------------------------------------------------
497	MParameterD pulpos("IntendedPulsePos");
498	// FIXME: Set a default which could be 1/3 of the digitization window
499	// pulpos.SetVal(40); // [ns]
500	plist.AddToList(&pulpos);
501
502	MParameterD trigger("TriggerPos");
503	trigger.SetVal(0);
504	plist.AddToList(&trigger);
505
506
507	// --------------------------------------------------------------------------------
508	// Setup container for the fix time offset, for residual time spread and for gapd time jitter
509	// --------------------------------------------------------------------------------
510	// Dominik and Sebastian on: fix time offsets
511	MMatrix fix_time_offsets_between_pixels_in_ns("MFixTimeOffset", 0);
512	plist.AddToList(&fix_time_offsets_between_pixels_in_ns);
513
514	// Jens Buss on: residual time spread
515	MParameterD resTimeSpread("ResidualTimeSpread");
516	resTimeSpread.SetVal(0.0);
517	plist.AddToList(&resTimeSpread);
518	// Jens Buss on: residual time spread
519	MParameterD gapdTimeJitter("GapdTimeJitter");
520	gapdTimeJitter.SetVal(0.0);
521	plist.AddToList(&gapdTimeJitter);
522
523	// --------------------------------------------------------------------------------
524	// Creation of binning objects and apply default settings
525	// --------------------------------------------------------------------------------
526	CreateBinningObjects(plist);
527
528	// --------------------------------------------------------------------------------
529	// --------------------------------------------------------------------------------
530	// Tasks Setup (Reading, Absorption, Reflector)
531	// --------------------------------------------------------------------------------
532	// --------------------------------------------------------------------------------
533
534	// --------------------------------------------------------------------------------
535	// Corsika read setup
536	// --------------------------------------------------------------------------------
537	MCorsikaRead read;
538	read.SetForceMode(fForceMode);
539
540	if (!seq.IsValid())
541	{
542	for (int i=0; i<args.GetNumArguments(); i++)
543	read.AddFile(args.GetArgumentStr(i));
544	}
545	else
546	read.AddFiles(iter);
547
548	if (read.GetNumFiles()==0 && fOperationMode==kModeData)
549	{
550	*fLog << err << "No files to read." << endl;
551	return kFALSE;
552	}
553
554	if (checkonly)
555	return read.CalcNumTotalEvents();
556
557	// --------------------------------------------------------------------------------
558	// Precut
559	// --------------------------------------------------------------------------------
560	MContinue precut("", "PreCut");
561	precut.IsInverted();
562	precut.SetAllowEmpty();
563
564	// --------------------------------------------------------------------------------
565	// MSimMMCS (don't know what this is doing) and calculation of the incident angle
566	// --------------------------------------------------------------------------------
567	MSimMMCS simmmcs;
568	const TString sin2 = "sin(MCorsikaEvtHeader.fZd)sin(MCorsikaRunHeader.fZdMinTMath::DegToRad())";
569	const TString cos2 = "cos(MCorsikaEvtHeader.fZd)cos(MCorsikaRunHeader.fZdMinTMath::DegToRad())";
570	const TString cos = "cos(MCorsikaEvtHeader.fAz-MCorsikaRunHeader.fAzMin*TMath::DegToRad())";
571
572	const TString form = "acos("+sin2+""+cos+"+"+cos2+")TMath::RadToDeg()";
573
574	MParameterCalc calcangle(form, "CalcIncidentAngle");
575	calcangle.SetNameParameter("IncidentAngle");
576
577	// --------------------------------------------------------------------------------
578	// Absorption tasks
579	// --------------------------------------------------------------------------------
580	// - Atmosphere (simatm)
581	// - photon detection efficiency (absapd)
582	// - mirror reflectivity (absmir)
583	// - angular acceptance of winston cones (cones)
584	// - transmission of winston cones (cones2)
585	// - additional photon acceptance (only motivated, not measured) (additionalPhotonAcceptance)
586	MSimAtmosphere simatm;
587	MSimAbsorption absapd("SimPhotonDetectionEfficiency");
588	MSimAbsorption absmir("SimMirrorReflectivity");
589	MSimAbsorption cones("SimConesAngularAcceptance");
590	MSimAbsorption cones2("SimConesTransmission");
591	MSimAbsorption additionalPhotonAcceptance("SimAdditionalPhotonAcceptance");
592	absapd.SetParName("PhotonDetectionEfficiency");
593	absmir.SetParName("MirrorReflectivity");
594	cones.SetParName("ConesAngularAcceptance");
595	cones.SetUseTheta();
596	cones2.SetParName("ConesTransmission");
597	additionalPhotonAcceptance.SetParName("AdditionalPhotonAcceptance");
598	additionalPhotonAcceptance.SetForce(kTRUE);
599
600	// --------------------------------------------------------------------------------
601	// Simulating pointing position and simulating reflector
602	// --------------------------------------------------------------------------------
603	MSimPointingPos pointing;
604	MSrcPosCalc srcposcam;
605
606	MSimReflector reflect;
607	reflect.SetNameGeomCam("GeomCones");
608	reflect.SetNameReflector("Reflector");
609	// MSimStarField stars;
610
611	// --------------------------------------------------------------------------------
612	// Continue tasks
613	// --------------------------------------------------------------------------------
614	// - Processing of the current events stops, if there aren't at least two photons (see cont3)
615	MContinue cont1("MPhotonEvent.GetNumPhotons<1", "ContEmpty1");
616	MContinue cont2("MPhotonEvent.GetNumPhotons<1", "ContEmpty2");
617	MContinue cont3("MPhotonEvent.GetNumPhotons<2", "ContEmpty3");
618	cont1.SetAllowEmpty();
619	cont2.SetAllowEmpty();
620	cont3.SetAllowEmpty();
621
622	// --------------------------------------------------------------------------------
623	// --------------------------------------------------------------------------------
624	// Tasks Setup (Cones, SiPM-Simulation, RandomPhotons, Camera, Trigger, DAQ)
625	// --------------------------------------------------------------------------------
626	// --------------------------------------------------------------------------------
627
628	// --------------------------------------------------------------------------------
629	// GeomApply, SimPSF, SimGeomCam, SimCalibrationSignal
630	// --------------------------------------------------------------------------------
631	// - MGeomApply only resizes all MGeomCam parameter container to the actual
632	// number of pixels
633	// - MSimPSF applies an additional smearing of the photons on the camera window
634	// on default it is switched of (set to 0), but can be applied by setting:
635	// MSimPSF.fSigma: <value>
636	// in the config file.
637	// Be aware, that there is also a smearing of the photons implemented in
638	// MSimReflector by setting:
639	// MReflector.SetSigmaPSF: <value>
640	// in the config file. This is at the moment done in the default config file.
641	// - MSimGeomCam identifies which pixel was hit by which photon
642	// - MSimCalibrationSignal is only used for simulated calibration runs
643	MGeomApply apply;
644
645	MSimPSF simpsf;
646
647	MSimGeomCam simgeom;
648	simgeom.SetNameGeomCam("GeomCones");
649
650	MSimCalibrationSignal simcal;
651	simcal.SetNameGeomCam("GeomCones");
652
653	// --------------------------------------------------------------------------------
654	// Simulation of random photons
655	// --------------------------------------------------------------------------------
656	// - be aware that here default values for the nsb rate and for the dark
657	// count rate are set. They will be overwritten if the values are defined
658	// in the config file
659	// - if you want to simulate a specific nsb rate you have to set the filename
660	// in the config file to an empty string and then you can specify the NSB rate:
661	// MSimRandomPhotons.FileNameNSB:
662	// MSimRandomPhotons.FrequencyNSB: 0.0
663
664	// Sky Quality Meter: 21.82M = 2.02e-4cd/m^2
665	// 1cd = 12.566 lm / 4pi sr
666	// FIXME: Simulate photons before cones and QE!
667	MSimRandomPhotons simnsb;
668	simnsb.SetFreq(5.8, 0.004); // ph/m^2/nm/sr/ns NSB, 4MHz dark count rate
669	simnsb.SetNameGeomCam("GeomCones");
670	// FIXME: How do we handle star-light? We need the rate but we also
671	// need to process it through the mirror!
672
673	// --------------------------------------------------------------------------------
674	// Simulation from the SiPM to the DAQ
675	// --------------------------------------------------------------------------------
676	// - MSimAPD simulates the whole behaviour of the SiPMs:
677	// (dead time of cells, crosstalk, afterpulses)
678	// - MSimExcessNoise adds a spread on the weight of each signal in the SiPMs
679	// - MSimBundlePhotons restructured the photon list of MPhotonEvent via a look
680	// up table. At the moment the tasks is skipped, due to the fact that there
681	// is no file name specified in the config file
682	// - MSimSignalCam only fills a SignalCam container for the cherenkov photons
683	// - MSimCamera simulates the behaviour of the camera:
684	// (electronic noise, accoupling, time smearing and offset for the photons,
685	// pulses injected by all photons)
686	// - MSimTrigger simulates the behaviour of the Trigger:
687	// (Adding patch voltages, applying clipping, applying discriminator, applying
688	// time over threshold and a possible trigger logic)
689	// - MContinue conttrig stops the processing of the current event if the trigger
690	// position is negativ (and therefore not valid)
691	// - MSimReadout simulates the behaviour of the readout:
692	// (Digitization and saturation of the ADC)
693	MSimAPD simapd;
694	simapd.SetNameGeomCam("GeomCones");
695
696	MSimExcessNoise simexcnoise;
697	MSimBundlePhotons simsum;
698	MSimSignalCam simsignal;
699	MSimCamera simcam;
700	MSimTrigger simtrig;
701	MContinue conttrig("TriggerPos.fVal<0", "ContTrigger");
702	MSimReadout simdaq;
703
704	// --------------------------------------------------------------------------------
705	// Standard analysis with hillas parameters for the true cherenkov photons
706	// --------------------------------------------------------------------------------
707	// Remove isolated pixels
708	MImgCleanStd clean(0, 0);
709	//clean.SetCleanLvl0(0); // The level above which isolated pixels are kept
710	clean.SetCleanRings(0);
711	clean.SetMethod(MImgCleanStd::kAbsolute);
712
713	//clean.SetNamePedPhotCam("MPedPhotFromExtractorRndm");
714
715	MHillasCalc hcalc;
716	hcalc.Disable(MHillasCalc::kCalcConc);
717
718
719	// --------------------------------------------------------------------------------
720	// Setup of histogram containers and the corresponding fill tasks
721	// --------------------------------------------------------------------------------
722	// Here is no functionality for the simulation, it is only visualization via
723	// showplot
724	MHn mhn1, mhn2, mhn3, mhn4;
725	SetupHist(mhn1);
726	SetupHist(mhn2);
727	SetupHist(mhn3);
728	SetupHist(mhn4);
729
730	MH3 mhtp("TriggerPos.fVal-IntendedPulsePos.fVal-PulseShape.GetWidth*1000/MRawRunHeader.GetFreqSampling");
731	mhtp.SetName("TrigPos");
732	mhtp.SetTitle("Trigger position w.r.t. the first photon (=Pulse.Xmin) hitting a detector [ns];#Delta T [ns]");
733
734	MH3 mhew("MPhotonStatistics.fLength");
735	mhew.SetName("TotLength");
736	mhew.SetTitle("Time between first and last photon hitting a detector;#Delta T [ns]");
737
738	MH3 mhed("MPhotonStatistics.fTimeMedDev");
739	mhed.SetName("MedLength");
740	mhed.SetTitle("Median deviation (1#sigma);#sigma(#Delta T) [ns]");
741
742	MFillH fillh1(&mhn1, "", "FillCorsika");
743	MFillH fillh2(&mhn2, "", "FillH2");
744	MFillH fillh3(&mhn3, "", "FillH3");
745	MFillH fillh4(&mhn4, "", "FillH4");
746	MFillH filltp(&mhtp, "", "FillTriggerPos");
747	MFillH fillew(&mhew, "", "FillEvtWidth");
748	MFillH filled(&mhed, "", "FillMedDev");
749	fillh1.SetNameTab("Corsika", "Distribution as simulated by Corsika");
750	fillh2.SetNameTab("Detectable", "Distribution of events hit the detector");
751	fillh3.SetNameTab("Triggered", "Distribution of triggered events");
752	fillh4.SetNameTab("Cleaned", "Distribution after cleaning and cuts");
753	filltp.SetNameTab("TrigPos", "TriggerPosition w.r.t the first photon");
754	fillew.SetNameTab("EvtWidth", "Time between first and last photon hitting a detector");
755	filled.SetNameTab("MedDev", "Median deviation of arrival time of photons hitting a detector");
756
757	MHPhotonEvent planeG(1, "HPhotonEventGround"); // Get from MaxImpact
758	MHPhotonEvent plane0(2, "HMirrorPlane0"); // Get from MReflector
759	//MHPhotonEvent plane1(2, "HMirrorPlane1");
760	MHPhotonEvent plane2(2, "HMirrorPlane2");
761	MHPhotonEvent plane3(2, "HMirrorPlane3");
762	MHPhotonEvent plane4(2, "HMirrorPlane4");
763	MHPhotonEvent planeF1(6, "HPhotonEventCamera"); // Get from MGeomCam
764	MHPhotonEvent planeF2(header.IsPointRun()?4:6, "HPhotonEventCones"); // Get from MGeomCam
765
766	plist.AddToList(&planeG);
767	plist.AddToList(&plane0);
768	plist.AddToList(&plane2);
769	plist.AddToList(&plane3);
770	plist.AddToList(&plane4);
771	plist.AddToList(&planeF1);
772	plist.AddToList(&planeF2);
773
774	//MHPSF psf;
775
776	MFillH fillG(&planeG, "MPhotonEvent", "FillGround");
777	MFillH fill0(&plane0, "MirrorPlane0", "FillReflector");
778	//MFillH fill1(&plane1, "MirrorPlane1", "FillCamShadow");
779	MFillH fill2(&plane2, "MirrorPlane2", "FillCandidates");
780	MFillH fill3(&plane3, "MirrorPlane3", "FillReflected");
781	MFillH fill4(&plane4, "MirrorPlane4", "FillFocal");
782	MFillH fillF1(&planeF1, "MPhotonEvent", "FillCamera");
783	MFillH fillF2(&planeF2, "MPhotonEvent", "FillCones");
784	//MFillH fillP(&psf, "MPhotonEvent", "FillPSF");
785
786	fillG.SetNameTab("Ground", "Photon distribution at ground");
787	fill0.SetNameTab("Reflector", "Photon distribution at reflector plane w/o camera shadow");
788	//fill1.SetNameTab("CamShadow", "Photon distribution at reflector plane w/ camera shadow");
789	fill2.SetNameTab("Candidates", "Can hit photon distribution at reflector plane w/ camera shadow");
790	fill3.SetNameTab("Reflected", "Photon distribution after reflector projected to reflector plane");
791	fill4.SetNameTab("Focal", "Photon distribution at focal plane");
792	fillF1.SetNameTab("Camera", "Photon distribution which hit the detector");
793	fillF2.SetNameTab("Cones", "Photon distribution after cones");
794	//fillP.SetNameTab("PSF", "Photon distribution after cones for all mirrors");
795
796
797	MHCamEvent evt0a(/10/3, "Signal", "Average signal (absolute);;S [ph]");
798	MHCamEvent evt0c(/10/3, "MaxSignal", "Maximum signal (absolute);;S [ph]");
799	MHCamEvent evt0d(/11/8, "ArrTm", "Time after first photon;;T [ns]");
800	evt0a.SetErrorSpread(kFALSE);
801	evt0c.SetCollectMax();
802
803	MContinue cut("", "Cut");
804
805	MFillH fillx0a(&evt0a, "MSignalCam", "FillAvgSignal");
806	MFillH fillx0c(&evt0c, "MSignalCam", "FillMaxSignal");
807	MFillH fillx0d(&evt0d, "MSignalCam", "FillArrTm");
808	MFillH fillx1("MHHillas", "MHillas", "FillHillas");
809	MFillH fillx3("MHHillasSrc", "MHillasSrc", "FillHillasSrc");
810	MFillH fillx4("MHNewImagePar", "MNewImagePar", "FillNewImagePar");
811	MFillH fillth("MHThreshold", "", "FillThreshold");
812	MFillH fillca("MHCollectionArea", "", "FillTrigArea");
813
814	fillth.SetNameTab("Threshold");
815	fillca.SetNameTab("TrigArea");
816
817	// --------------------------------------------------------------------------------
818	// Formating of the outputfilepath
819	// --------------------------------------------------------------------------------
820	if (!fFileOut.IsNull())
821	{
822	const Ssiz_t dot = fFileOut.Last('.');
823	const Ssiz_t slash = fFileOut.Last('/');
824	if (dot>slash)
825	fFileOut = fFileOut.Remove(dot);
826	}
827
828	const char *regex = fRunNumber<0 ?
829	"s/[cC][eE][rR]([0-9]+)(-[tT][eE][lL]([0-9]+))?/%s\\/00$1.%03d%%s_MonteCarlo$2.root/" :
830	"s/[cC][eE][rR][0-9]+([0-9][0-9][0-9])(-[tT][eE][lL]([0-9]+))?/%s\\/%08d.$1%%s_MonteCarlo$2.root/";
831
832	const char *fmt = fFileOut.IsNull() ?
833	Form(regex, Esc(fPathOut).Data(), header.GetRunNumber()) :
834	Form("%s/%s%%s.root", Esc(fPathOut).Data(), Esc(fFileOut).Data());
835
836	const TString rule1(Form(fmt, fFileOut.IsNull()?"_R":""));
837	const TString rule2(Form(fmt, "_Y"));
838	const TString rule4(Form(fmt, "_I"));
839	TString rule3(Form(fmt, Form("_%c", header.GetRunTypeChar())));
840
841	// --------------------------------------------------------------------------------
842	// Setup of the outputfile tasks
843	// --------------------------------------------------------------------------------
844	MWriteRootFile write4a( 2, rule4, fOverwrite?"RECREATE":"NEW", "Star file");
845	MWriteRootFile write4b( 2, rule4, fOverwrite?"RECREATE":"NEW", "Star file");
846	MWriteRootFile write3b( 2, rule3, fOverwrite?"RECREATE":"NEW", "Camera file");
847	MWriteRootFile write2a( 2, rule2, fOverwrite?"RECREATE":"NEW", "Signal file");
848	MWriteRootFile write2b( 2, rule2, fOverwrite?"RECREATE":"NEW", "Signal file");
849	MWriteRootFile write1a( 2, rule1, fOverwrite?"RECREATE":"NEW", "Reflector file");
850	MWriteRootFile write1b( 2, rule1, fOverwrite?"RECREATE":"NEW", "Reflector file");
851
852	if (fWriteFitsFile)
853	rule3.ReplaceAll(".root", ".fits");
854
855	MWriteFitsFile write3af( 2, rule3, fOverwrite?"RECREATE":"NEW", "Camera file");
856
857	if (fWriteFitsFile)
858	rule3 = "/dev/null";
859
860	MWriteRootFile write3ar( 2, rule3, fOverwrite?"RECREATE":"NEW", "Camera file");
861
862	MTask &write3a = fWriteFitsFile ? static_cast<MTask&>(write3af) : static_cast<MTask&>(write3ar);
863
864	//SetupHeaderKeys(write3af,header);
865	SetupVetoColumns(write3af);
866
867	write3af.SetBytesPerSample("Data", 2);
868
869	write1a.SetName("WriteRefData");
870	write1b.SetName("WriteRefMC");
871	write2a.SetName("WriteSigData");
872	write2b.SetName("WriteSigMC");
873	write3a.SetName("WriteCamData");
874	write3b.SetName("WriteCamMC");
875	write4a.SetName("WriteImgData");
876	write4b.SetName("WriteImgMC");
877
878	SetupCommonFileStructure(write1a);
879	SetupCommonFileStructure(write2a);
880	SetupCommonFileStructure(write3ar);
881	SetupCommonFileStructure(write3af);
882	SetupCommonFileStructure(write4a);
883
884	// R: Dedicated file structure
885	write1a.AddContainer("MPhotonEvent", "Events");
886
887	// Y: Dedicated file structure
888	write2a.AddContainer("MPedPhotFromExtractorRndm", "RunHeaders", kTRUE, 1); // FIXME: Needed for the signal files to be display in MARS
889	write2a.AddContainer("MSignalCam", "Events");
890
891	// D: Dedicated file structure
892	write3af.AddContainer("ElectronicNoise", "RunHeaders", kTRUE, 1);
893	write3af.AddContainer("IntendedPulsePos", "RunHeaders", kTRUE, 1);
894	write3af.AddContainer("ResidualTimeSpread", "RunHeaders", kTRUE, 1);
895	write3af.AddContainer("GapdTimeJitter", "RunHeaders", kTRUE, 1);
896	write3af.AddContainer("MRawEvtData", "Events");
897	write3af.AddContainer("MTruePhotonsPerPixelCont", "Events");
898
899	write3ar.AddContainer("ElectronicNoise", "RunHeaders", kTRUE, 1);
900	write3ar.AddContainer("IntendedPulsePos", "RunHeaders", kTRUE, 1);
901	write3ar.AddContainer("MRawEvtData", "Events");
902	// It doesn't make much sene to write this information
903	// to the file because the units are internal units not
904	// related to the output samples
905	// if (header.IsDataRun() \|\| fForceTrigger)
906	// write3a.AddContainer("TriggerPos", "Events");
907
908	// I: Dedicated file structure
909	write4a.AddContainer("MHillas", "Events");
910	write4a.AddContainer("MHillasSrc", "Events");
911	write4a.AddContainer("MImagePar", "Events");
912	write4a.AddContainer("MNewImagePar", "Events");
913
914	// Basic MC data
915	write1b.AddContainer("MMcEvtBasic", "OriginalMC");
916	write2b.AddContainer("MMcEvtBasic", "OriginalMC");
917	write3b.AddContainer("MMcEvtBasic", "OriginalMC");
918	write4b.AddContainer("MMcEvtBasic", "OriginalMC");
919
920	// --------------------------------------------------------------------------------
921	// --------------------------------------------------------------------------------
922	// Filling of tasks in tasklist
923	// --------------------------------------------------------------------------------
924	// --------------------------------------------------------------------------------
925
926	if (header.IsDataRun())
927	{
928	tasks.AddToList(&read); // Reading Corsika
929	tasks.AddToList(&precut); // Precut
930	tasks.AddToList(&pointing); // Simulating pointing
931	tasks.AddToList(&srcposcam); // calculate origin in camera
932	tasks.AddToList(&simmmcs); // Simulating MMCS
933	if (!fPathOut.IsNull() && !HasNullOut()) // Write Tasks for corsika infos
934	{
935	//tasks.AddToList(&write1b);
936	tasks.AddToList(&write2b);
937	if (fCamera)
938	tasks.AddToList(&write3b);
939	if (header.IsDataRun())
940	tasks.AddToList(&write4b);
941	}
942	// if (header.IsPointRun())
943	// tasks.AddToList(&stars);
944	if (1)
945	tasks.AddToList(&simatm); // Here because before fillh1 // atmosphere absorption
946	tasks.AddToList(&calcangle); // calculation incident angle
947	tasks.AddToList(&fillh1); // fill histogram task
948	tasks.AddToList(&fillG); // fill histogram task
949	if (!header.IsPointRun())
950	{
951	tasks.AddToList(&absapd); // photon detection efficiency of the apds
952	tasks.AddToList(&absmir); // mirror reflectivity
953	tasks.AddToList(&additionalPhotonAcceptance); // addition photon acceptance
954	if (0)
955	tasks.AddToList(&simatm); // FASTER? // be aware this is 'commented'
956	}
957	tasks.AddToList(&reflect); // Simulation of the reflector
958	if (!header.IsPointRun())
959	{
960	tasks.AddToList(&fill0); // fill histogram task
961	//tasks.AddToList(&fill1);
962	tasks.AddToList(&fill2); // fill histogram task
963	tasks.AddToList(&fill3); // fill histogram task
964	tasks.AddToList(&fill4); // fill histogram task
965	tasks.AddToList(&fillF1); // fill histogram task
966	}
967	tasks.AddToList(&cones); // angular acceptance of winston cones
968	tasks.AddToList(&cones2); // transmission of winston cones
969	//if (header.IsPointRun())
970	// tasks.AddToList(&fillP);
971	tasks.AddToList(&cont1); // continue if at least 2 photons
972	}
973	// -------------------------------
974	tasks.AddToList(&apply); // apply geometry to MGeomCam containers
975	if (header.IsDataRun())
976	{
977	tasks.AddToList(&simpsf); // simulates additional psf (NOT used in default mode)
978	tasks.AddToList(&simgeom); // tag which pixel is hit by which photon
979	tasks.AddToList(&cont2); // continue if at least 2 photons
980	if (!header.IsPointRun())
981	{
982	tasks.AddToList(&fillF2); // fill histogram task
983	tasks.AddToList(&fillh2); // fill histogram task
984	}
985	tasks.AddToList(&cont3); // continue if at least 3 photons
986	}
987	if (fCamera)
988	{
989	if (header.IsPedestalRun() \|\| header.IsCalibrationRun())
990	tasks.AddToList(&simcal); // add calibration signal for calibration runs
991	tasks.AddToList(&simnsb); // simulate nsb and dark counts
992	tasks.AddToList(&simapd); // simulate SiPM behaviour (dead time, crosstalk ...)
993	tasks.AddToList(&simexcnoise); // add excess noise
994	}
995	tasks.AddToList(&simsum); // bundle photons (NOT used in default mode)
996	if (fCamera)
997	{
998	tasks.AddToList(&simcam); // simulate camera behaviour (creates analog signal)
999	if (header.IsDataRun() \|\| fForceTrigger)
1000	tasks.AddToList(&simtrig); // simulate trigger
1001	tasks.AddToList(&conttrig); // continue if trigger pos is valid
1002	tasks.AddToList(&simdaq); // simulate data aquisition
1003	}
1004	tasks.AddToList(&simsignal); // What do we do if signal<0? // fill MSimSignal container
1005	if (!fPathOut.IsNull() && !HasNullOut())
1006	{
1007	//tasks.AddToList(&write1a);
1008	if (!header.IsPedestalRun())
1009	tasks.AddToList(&write2a); // outputtask
1010	if (fCamera)
1011	tasks.AddToList(&write3a); // outputtask (this is the raw data writing tasks)
1012	}
1013	// -------------------------------
1014	if (fCamera)
1015	{
1016	// FIXME: MHCollectionArea Trigger Area!
1017	if (header.IsDataRun())
1018	tasks.AddToList(&fillh3); // fill histogram task
1019	tasks.AddToList(&filltp); // fill histogram task
1020	}
1021	if (header.IsDataRun())
1022	{
1023	tasks.AddToList(&fillew); // fill histogram task
1024	tasks.AddToList(&filled); // fill histogram task
1025	}
1026	if (!header.IsPedestalRun())
1027	{
1028	tasks.AddToList(&fillx0a); // fill histogram task
1029	tasks.AddToList(&fillx0c); // fill histogram task
1030	}
1031	if (header.IsDataRun())
1032	{
1033	tasks.AddToList(&clean); // Cleaning for standard analysis
1034	tasks.AddToList(&hcalc); // Hillas parameter calculation
1035	tasks.AddToList(&cut);
1036	tasks.AddToList(&fillx0d); // fill histogram task
1037	tasks.AddToList(&fillx1); // fill histogram task
1038	//tasks.AddToList(&fillx2);
1039	tasks.AddToList(&fillx3); // fill histogram task
1040	tasks.AddToList(&fillx4); // fill histogram task
1041	if (!HasNullOut())
1042	tasks.AddToList(&write4a);
1043	//tasks.AddToList(&fillx5);
1044
1045	tasks.AddToList(&fillh4); // fill histogram task
1046	tasks.AddToList(&fillth); // fill histogram task
1047	tasks.AddToList(&fillca); // fill histogram task
1048	}
1049
1050
1051	// --------------------------------------------------------------------------------
1052	// --------------------------------------------------------------------------------
1053	// Event loop and processing display
1054	// --------------------------------------------------------------------------------
1055	// --------------------------------------------------------------------------------
1056
1057	// --------------------------------------------------------------------------------
1058	// Creation of event loop
1059	// --------------------------------------------------------------------------------
1060	tasks.SetAccelerator(MTask::kAccDontReset\|MTask::kAccDontTime);
1061
1062	// Create and setup the eventloop
1063	MEvtLoop evtloop(fName);
1064	evtloop.SetParList(&plist);
1065	evtloop.SetDisplay(fDisplay);
1066	evtloop.SetLogStream(&gLog);
1067	if (!SetupEnv(evtloop))
1068	return kFALSE;
1069
1070	// FIXME: If pedestal file given use the values from this file
1071
1072
1073	// --------------------------------------------------------------------------------
1074	// Preparation of the graphicalk display which is shown while processing
1075	// --------------------------------------------------------------------------------
1076
1077	MGeomCam cam = static_cast<MGeomCam>(env2.GetCont());
1078
1079	write3af.SetGeometry(cam->ClassName());
1080
1081	MBinning binsd = (MBinning) plist.FindObject("BinningDistC");
1082	MBinning binsd0 = (MBinning) plist.FindObject("BinningDist");
1083	MBinning binstr = (MBinning) plist.FindObject("BinningTrigPos");
1084	if (binstr->IsDefault())
1085	binstr->SetEdgesLin(150, -header.GetFreqSampling()/1000.*shape.GetWidth(),
1086	header.GetFreqSampling()/1000.*(header.GetNumSamples()+shape.GetWidth()));
1087
1088	if (binsd->IsDefault() && cam)
1089	binsd->SetEdgesLin(100, 0, cam->GetMaxRadius()*cam->GetConvMm2Deg());
1090
1091	if (binsd0->IsDefault() && cam)
1092	binsd0->SetEdgesLin(100, 0, cam->GetMaxRadius()*cam->GetConvMm2Deg());
1093
1094
1095	header.Print();
1096
1097	// FIXME: Display acceptance curves!
1098
1099	if (fDisplay)
1100	{
1101	TCanvas &c = fDisplay->AddTab("Info");
1102	c.Divide(2,2);
1103
1104	c.cd(1);
1105	gPad->SetBorderMode(0);
1106	gPad->SetFrameBorderMode(0);
1107	gPad->SetGridx();
1108	gPad->SetGridy();
1109	gROOT->SetSelectedPad(0);
1110	shape.DrawClone()->SetBit(kCanDelete);
1111
1112	if (env0.GetCont() && (header.IsDataRun() \|\| header.IsPointRun()))
1113	{
1114	c.cd(3);
1115	gPad->SetBorderMode(0);
1116	gPad->SetFrameBorderMode(0);
1117	gPad->SetGridx();
1118	gPad->SetGridy();
1119	gROOT->SetSelectedPad(0);
1120	static_cast<MOptics*>(env0.GetCont())->DrawClone("line")->SetBit(kCanDelete);
1121	}
1122
1123	if (fCamera)
1124	{
1125	if (env1.GetCont())
1126	{
1127	c.cd(2);
1128	gPad->SetBorderMode(0);
1129	gPad->SetFrameBorderMode(0);
1130	gPad->SetGridx();
1131	gPad->SetGridy();
1132	gROOT->SetSelectedPad(0);
1133	MHCamera hcam = new MHCamera(static_cast<MGeomCam&>(env1.GetCont()));
1134	hcam->SetStats(kFALSE);
1135	hcam->SetBit(MHCamera::kNoLegend);
1136	hcam->SetBit(kCanDelete);
1137	hcam->Draw();
1138	}
1139
1140	if (cam)
1141	{
1142	c.cd(4);
1143	gPad->SetBorderMode(0);
1144	gPad->SetFrameBorderMode(0);
1145	gPad->SetGridx();
1146	gPad->SetGridy();
1147	gROOT->SetSelectedPad(0);
1148	MHCamera hcam = new MHCamera(cam);
1149	hcam->SetStats(kFALSE);
1150	hcam->SetBit(MHCamera::kNoLegend);
1151	hcam->SetBit(kCanDelete);
1152	hcam->Draw();
1153	}
1154	}
1155
1156	TCanvas &d = fDisplay->AddTab("Info2");
1157	d.Divide(2,3);
1158
1159	d.cd(1);
1160	gPad->SetBorderMode(0);
1161	gPad->SetFrameBorderMode(0);
1162	gPad->SetGridx();
1163	gPad->SetGridy();
1164	gROOT->SetSelectedPad(0);
1165	splinepde.DrawClone()->SetBit(kCanDelete);
1166
1167	d.cd(3);
1168	gPad->SetBorderMode(0);
1169	gPad->SetFrameBorderMode(0);
1170	gPad->SetGridx();
1171	gPad->SetGridy();
1172	gROOT->SetSelectedPad(0);
1173	splinecones2.DrawClone()->SetBit(kCanDelete);
1174
1175	d.cd(5);
1176	gPad->SetBorderMode(0);
1177	gPad->SetFrameBorderMode(0);
1178	gPad->SetGridx();
1179	gPad->SetGridy();
1180	gROOT->SetSelectedPad(0);
1181	splinemirror.DrawClone()->SetBit(kCanDelete);
1182
1183	d.cd(2);
1184	gPad->SetBorderMode(0);
1185	gPad->SetFrameBorderMode(0);
1186	gPad->SetGridx();
1187	gPad->SetGridy();
1188	gROOT->SetSelectedPad(0);
1189	splinecones.DrawClone()->SetBit(kCanDelete);
1190	//splinecones2.DrawClone("same")->SetBit(kCanDelete);
1191
1192	d.cd(6);
1193	gPad->SetBorderMode(0);
1194	gPad->SetFrameBorderMode(0);
1195	gPad->SetGridx();
1196	gPad->SetGridy();
1197	gROOT->SetSelectedPad(0);
1198	MParSpline sall = (MParSpline)splinepde.DrawClone();
1199	//all->SetTitle("Combined acceptance");
1200	sall->SetBit(kCanDelete);
1201	if (splinemirror.GetSpline())
1202	sall->Multiply(*splinemirror.GetSpline());
1203	if (splinecones2.GetSpline())
1204	sall->Multiply(*splinecones2.GetSpline());
1205	}
1206
1207	// --------------------------------------------------------------------------------
1208	// Perform event loop
1209	// --------------------------------------------------------------------------------
1210
1211	// Execute first analysis
1212	if (!evtloop.Eventloop(fMaxEvents))
1213	{
1214	gLog << err << GetDescriptor() << ": Failed." << endl;
1215	return kFALSE;
1216	}
1217
1218	//-------------------------------------------
1219	// FIXME: Display Spline in tab
1220	// FIXME: Display Reflector in tab
1221	// FIXME: Display Routing(?) in tab
1222	// FIXME: Display Camera(s) in tab
1223	//-------------------------------------------
1224
1225	if (!WriteResult(plist, seq, header.GetRunNumber()))
1226	return kFALSE;
1227
1228	*fLog << all << GetDescriptor() << ": Done." << endl << endl << endl;;
1229
1230	return kTRUE;
1231	}

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

Context Navigation

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

Download in other formats: