Context Navigation

source: branches/Corsika7500Compatibility/mjobs/MJSimulation.cc@ 18454

Visit:

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

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

Download in other formats: