Context Navigation

source: trunk/MagicSoft/Mars/mjobs/MJSimulation.cc@ 9439

Visit:

Last change on this file since 9439 was 9428, checked in by tbretz, 16 years ago
* empty log message *
File size: 27.9 KB

Line
1	/* ======================================================================== *\
2	!
3	! *
4	! * This file is part of MARS, the MAGIC Analysis and Reconstruction
5	! * Software. It is distributed to you in the hope that it can be a useful
6	! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
7	! * It is distributed WITHOUT ANY WARRANTY.
8	! *
9	! * Permission to use, copy, modify and distribute this software and its
10	! * documentation for any purpose is hereby granted without fee,
11	! * provided that the above copyright notice appear in all copies and
12	! * that both that copyright notice and this permission notice appear
13	! * in supporting documentation. It is provided "as is" without express
14	! * or implied warranty.
15	! *
16	!
17	!
18	! Author(s): Thomas Bretz, 1/2009 <mailto:tbretz@astro.uni-wuerzburg.de>
19	!
20	! Copyright: MAGIC Software Development, 2000-2009
21	!
22	!
23	\* ======================================================================== */
24
25	/////////////////////////////////////////////////////////////////////////////
26	//
27	// MJSimulation
28	//
29	//
30	// Force reading a corsika file even if the footer (RUNE-section) is missing
31	// by setting fForceMode to kTRUE or from the resource file by
32	//
33	// ForceMode: Yes
34	//
35	//
36	// To switch off the simulation of the camera electronics, use:
37	//
38	// Camera: Off
39	//
40	// Note, that the border of the camera and the propertied of the cones
41	// are still simulated (simply because it is fast). Furthermore, this
42	// switches off the trigger for the output, i.e. all data which deposits
43	// at least one photon in at least one pixel is written to the output file.
44	//
45	//
46	// In case of a pedestal or calibration run the artificial trigger can
47	// be "switched off" and the cosmics trrigger "switched on" by setting
48	// fForceTrigger to kTRUE or from the resource file by
49	//
50	// ForceTrigger: Yes
51	//
52	//
53	/////////////////////////////////////////////////////////////////////////////
54	#include "MJSimulation.h"
55
56	#include <TEnv.h>
57	#include <TCanvas.h>
58
59	// Core
60	#include "MLog.h"
61	#include "MLogManip.h"
62
63	#include "MArgs.h"
64	//#include "MDirIter.h"
65	#include "MParList.h"
66	#include "MTaskList.h"
67	#include "MEvtLoop.h"
68
69	#include "MStatusDisplay.h"
70
71	// Tasks
72	#include "MCorsikaRead.h"
73	#include "MContinue.h"
74	#include "MFillH.h"
75	#include "MGeomApply.h"
76	#include "MHillasCalc.h"
77	#include "MImgCleanStd.h"
78	#include "MWriteRootFile.h"
79
80	#include "MSimMMCS.h"
81	#include "MSimAbsorption.h"
82	#include "MSimAtmosphere.h"
83	#include "MSimReflector.h"
84	#include "MSimPointingPos.h"
85	#include "MSimPSF.h"
86	#include "MSimGeomCam.h"
87	#include "MSimSignalCam.h"
88	#include "MSimAPD.h"
89	#include "MSimExcessNoise.h"
90	#include "MSimCamera.h"
91	#include "MSimTrigger.h"
92	#include "MSimReadout.h"
93	#include "MSimRandomPhotons.h"
94	#include "MSimBundlePhotons.h"
95	#include "MSimCalibrationSignal.h"
96
97	// Histograms
98	#include "MBinning.h"
99
100	#include "MHn.h"
101	#include "MHCamera.h"
102	#include "MHCamEvent.h"
103	#include "MHPhotonEvent.h"
104
105	// Container
106	#include "MRawRunHeader.h"
107	#include "MParameters.h"
108	#include "MReflector.h"
109	#include "MParEnv.h"
110	#include "MSpline3.h"
111	#include "MParSpline.h"
112	#include "MGeomCam.h"
113
114	#include "MPedestalCam.h"
115	#include "MPedestalPix.h"
116
117	ClassImp(MJSimulation);
118
119	using namespace std;
120
121	// --------------------------------------------------------------------------
122	//
123	// Default constructor.
124	//
125	// Sets fRuns to 0, fExtractor to NULL, fDataCheck to kFALSE
126	//
127	MJSimulation::MJSimulation(const char name, const char title)
128	: fForceMode(kFALSE), fCamera(kTRUE) ,fForceTrigger(kFALSE)
129	{
130	fName = name ? name : "MJSimulation";
131	fTitle = title ? title : "Standard analysis and reconstruction";
132	}
133
134	Bool_t MJSimulation::CheckEnvLocal()
135	{
136	fForceMode = GetEnv("ForceMode", fForceMode);
137	fForceTrigger = GetEnv("ForceTrigger", fForceTrigger);
138	fCamera = GetEnv("Camera", fCamera);
139
140	return kTRUE;
141	}
142
143	Bool_t MJSimulation::WriteResult(const MParList &plist)
144	{
145	if (fPathOut.IsNull())
146	{
147	*fLog << inf << "No output path specified via SetPathOut - no output written." << endl;
148	return kTRUE;
149	}
150
151	TObjArray cont;
152	cont.Add(const_cast<TEnv*>(GetEnv()));
153	//cont.Add(const_cast<MSequence*>(&fSequence));
154
155	cont.Add(plist.FindObject("PulseShape"));
156	cont.Add(plist.FindObject("Reflector"));
157	cont.Add(plist.FindObject("MGeomCam"));
158	cont.Add(plist.FindObject("GeomCones"));
159
160	if (fDisplay)
161	{
162	// TString title = "-- Reflector: ";
163	// title += fSequence.GetSequence();
164	// title += " --";
165	// fDisplay->SetTitle(title, kFALSE);
166	fDisplay->SetTitle("Ceres", kFALSE);
167
168	cont.Add(fDisplay);
169	}
170
171	// const TString oname = Form("reflector%08d.root", fSequence.GetSequence());
172	const TString oname = "ceres.root";
173	return WriteContainer(cont, oname, "RECREATE");
174	}
175
176	void MJSimulation::SetupHist(MHn &hist) const
177	{
178	hist.AddHist("MCorsikaEvtHeader.fTotalEnergy");
179	hist.InitName("Energy");
180	hist.InitTitle("Energy;E [GeV]");
181	hist.SetLog(kTRUE, kTRUE, kFALSE);
182
183	hist.AddHist("MPhotonEvent.GetNumExternal");
184	hist.InitName("Size");
185	hist.InitTitle("Size;S [#]");
186	hist.SetLog(kTRUE, kTRUE, kFALSE);
187
188	/*
189	hist.AddHist("-MCorsikaEvtHeader.fX/100","-MCorsikaEvtHeader.fY/100");
190	hist.SetDrawOption("colz");
191	hist.InitName("Impact;Impact;Impact");
192	hist.InitTitle("Impact;West <--> East [m];South <--> North [m]");
193	hist.SetAutoRange();
194	*/
195
196	hist.AddHist("MCorsikaEvtHeader.GetImpact/100");
197	hist.InitName("Impact");
198	hist.InitTitle("Impact;Impact [m]");
199	hist.SetAutoRange();
200
201	hist.AddHist("MCorsikaEvtHeader.fFirstInteractionHeight/100000");
202	hist.InitName("Height");
203	hist.InitTitle("FirstInteractionHeight;h [km]");
204
205	hist.AddHist("(MCorsikaEvtHeader.fAz+MCorsikaRunHeader.fMagneticFieldAz)TMath::RadToDeg()", "MCorsikaEvtHeader.fZdTMath::RadToDeg()");
206	hist.InitName("SkyOrigin;Az;Zd");
207	hist.InitTitle("Sky Origin;Az [\\deg];Zd [\\deg]");
208	hist.SetDrawOption("colz");
209	hist.SetAutoRange();
210
211	TString sin2 = "sin(MCorsikaEvtHeader.fZd)sin(MCorsikaRunHeader.fZdMinTMath::DegToRad())";
212	TString cos2 = "cos(MCorsikaEvtHeader.fZd)cos(MCorsikaRunHeader.fZdMinTMath::DegToRad())";
213	TString cos = "cos(MCorsikaEvtHeader.fAz-MCorsikaRunHeader.fAzMin*TMath::DegToRad())";
214
215	TString form = "acos("+sin2+""+cos+"+"+cos2+")TMath::RadToDeg()";
216
217	hist.AddHist(form);
218	hist.InitName("ViewCone");
219	hist.InitTitle("Incident Angle;\\alpha [\\deg]");
220	}
221
222	void MJSimulation::SetupCommonFileStructure(MWriteRootFile &write) const
223	{
224	// Common run headers
225	write.AddContainer("MMcCorsikaRunHeader", "RunHeaders", kFALSE);
226	write.AddContainer("MCorsikaRunHeader", "RunHeaders", kFALSE);
227	write.AddContainer("MRawRunHeader", "RunHeaders");
228	write.AddContainer("MGeomCam", "RunHeaders");
229	write.AddContainer("MMcRunHeader", "RunHeaders");
230
231	// Common events
232	write.AddContainer("MCorsikaEvtHeader", "Events", kFALSE);
233	write.AddContainer("MRawEvtHeader", "Events");
234	write.AddContainer("MMcEvt", "Events");
235	}
236
237	Bool_t MJSimulation::Process(const MArgs &args)
238	{
239	/*
240	if (!fSequence.IsValid())
241	{
242	*fLog << err << "ERROR - Sequence invalid!" << endl;
243	return kFALSE;
244	}
245	*/
246
247	//if (!HasWritePermission(GetPathOut()))
248	// return kFALSE;
249
250	*fLog << inf;
251	fLog->Separator(GetDescriptor());
252
253	if (!CheckEnv())
254	return kFALSE;
255
256	*fLog << warn << "FIXME: Monte Carlo simulation: Sequences not supported yet.";
257	//*fLog << fSequence.GetFileName() << endl;
258	*fLog << endl;
259
260	// --------------------------------------------------------------------------------
261
262	//MDirIter iter;
263	//if (fSequence.GetRuns(iter, MSequence::kCalibrated)<=0)
264	// return kFALSE;
265
266	// Setup Parlist
267	MParList plist;
268	plist.AddToList(this); // take care of fDisplay!
269
270	// setup TaskList
271	MTaskList tasks;
272	plist.AddToList(&tasks);
273
274	// --------------------------------------------------------------------------------
275
276	// FIXME: Allow empty GeomCones!
277	MParEnv env0("Reflector");
278	MParEnv env1("GeomCones");
279	MParEnv env2("MGeomCam");
280	env0.SetClassName("MReflector");
281	env1.SetClassName("MGeomCam");
282	env2.SetClassName("MGeomCam");
283	plist.AddToList(&env0);
284	plist.AddToList(&env1);
285	plist.AddToList(&env2);
286
287	plist.FindCreateObj("MPedPhotCam", "MPedPhotFromExtractorRndm");
288
289	MParSpline shape("PulseShape");
290	plist.AddToList(&shape);
291
292	// * FIXME * FIXME * FIXME *
293	plist.FindCreateObj("MMcRunHeader");
294
295	MRawRunHeader header;
296	header.SetValidMagicNumber();
297	//header.InitFadcType(3);
298
299	header.SetRunType(MRawRunHeader::kRTMonteCarlo\|MRawRunHeader::kRTData);
300	if (args.GetNumArguments()==1)
301	{
302	if (!args.GetArgumentStr(0).CompareTo("pedestal", TString::kIgnoreCase))
303	{
304	header.SetRunType(MRawRunHeader::kRTMonteCarlo\|MRawRunHeader::kRTPedestal);
305	header.SetSourceInfo("Pedestal");
306	}
307	if (!args.GetArgumentStr(0).CompareTo("calibration", TString::kIgnoreCase))
308	{
309	header.SetRunType(MRawRunHeader::kRTMonteCarlo\|MRawRunHeader::kRTCalibration);
310	header.SetSourceInfo("Calibration");
311	}
312	if (!args.GetArgumentStr(0).CompareTo("pointrun", TString::kIgnoreCase))
313	header.SetRunType(MRawRunHeader::kRTMonteCarlo\|MRawRunHeader::kRTPointRun);
314	}
315
316	// FIXME: Move to MSimPointingPos, MSimCalibrationSignal
317	// Can we use this as input for MSimPointingPos?
318	header.SetObservation("On", "MonteCarlo");
319	plist.AddToList(&header);
320	// ++++++++ FIXME FIXME FIXME +++++++++++++
321
322	/*
323	MPedestalCam pedcam;
324	pedcam.Init(geomcam.GetNumPixels());
325	for (UInt_t i=0; i<geomcam.GetNumPixels(); i++)
326	pedcam[i].Set(128./header.GetScale(), 22.5/header.GetScale());
327	plist.AddToList(&pedcam);
328	*/
329
330	// -------------------------------------------------------------------
331
332	MCorsikaRead read;
333	read.SetForceMode(fForceMode);
334
335	for (int i=0; i<args.GetNumArguments(); i++)
336	read.AddFile(args.GetArgumentStr(i));
337
338	MSimMMCS simmmcs;
339
340	MParSpline splinepde("PhotonDetectionEfficiency");
341	MParSpline splinemirror("MirrorReflectivity");
342	MParSpline splinecones("ConesAngularAcceptance");
343	plist.AddToList(&splinepde);
344	plist.AddToList(&splinemirror);
345	plist.AddToList(&splinecones);
346
347	MSimAtmosphere simatm;
348	MSimAbsorption absapd("SimPhotonDetectionEfficiency");
349	MSimAbsorption absmir("SimMirrorReflectivity");
350	MSimAbsorption cones("SimConesAngularAcceptance");
351	absapd.SetParName("PhotonDetectionEfficiency");
352	absmir.SetParName("MirrorReflectivity");
353	cones.SetParName("ConesAngularAcceptance");
354	cones.SetUseTheta();
355
356	MSimPointingPos pointing;
357
358	MSimReflector reflect;
359	reflect.SetNameGeomCam("GeomCones");
360	// MSimStarField stars;
361
362	MContinue cont1("MPhotonEvent.GetNumPhotons<1", "ContEmpty1");
363	MContinue cont2("MPhotonEvent.GetNumPhotons<1", "ContEmpty2");
364	MContinue cont3("MPhotonEvent.GetNumPhotons<2", "ContEmpty3");
365
366	// -------------------------------------------------------------------
367
368	MBinning binse( 100, 1, 100000, "BinningEnergy", "log");
369	MBinning binsth( 70, 0.9, 900000, "BinningThreshold", "log");
370	MBinning binsee( 35, 0.9, 900000, "BinningEnergyEst", "log");
371	MBinning binss( 100, 1, 10000000, "BinningSize", "log");
372	// MBinning binsi( 100, -500, 500, "BinningImpact");
373	MBinning binsi( 55, 0, 1100, "BinningImpact");
374	MBinning binsh( 150, 0, 50, "BinningHeight");
375	MBinning binsaz(720, -360, 360, "BinningAz");
376	MBinning binszd( 70, 0, 70, "BinningZd");
377	MBinning binsvc( 45, 0, 9, "BinningViewCone");
378	MBinning binsel(150, 0, 50, "BinningTotLength");
379	MBinning binsew(150, 0, 15, "BinningMedLength");
380	MBinning binsd("BinningDistC");
381	MBinning binsd0("BinningDist");
382	MBinning binstr("BinningTrigPos");
383
384	plist.AddToList(&binsee);
385	plist.AddToList(&binse);
386	plist.AddToList(&binsth);
387	plist.AddToList(&binss);
388	plist.AddToList(&binsi);
389	plist.AddToList(&binsh);
390	plist.AddToList(&binszd);
391	plist.AddToList(&binsaz);
392	plist.AddToList(&binsvc);
393	plist.AddToList(&binsel);
394	plist.AddToList(&binsew);
395	plist.AddToList(&binstr);
396	plist.AddToList(&binsd);
397	plist.AddToList(&binsd0);
398
399	MHn mhn1, mhn2, mhn3, mhn4;
400	SetupHist(mhn1);
401	SetupHist(mhn2);
402	SetupHist(mhn3);
403	SetupHist(mhn4);
404
405	MH3 mhtp("TriggerPos.fVal-IntendedPulsePos.fVal-PulseShape.GetWidth");
406	mhtp.SetName("TrigPos");
407	mhtp.SetTitle("Trigger position w.r.t. the first photon hitting a detector");
408
409	MH3 mhew("MPhotonStatistics.fLength");
410	mhew.SetName("TotLength");
411	mhew.SetTitle("Time between first and last photon hitting a detector;L [ns]");
412
413	MH3 mhed("MPhotonStatistics.fTimeMedDev");
414	mhed.SetName("MedLength");
415	mhed.SetTitle("Median deviation (1\\sigma);L [ns]");
416
417	MFillH fillh1(&mhn1, "", "FillCorsika");
418	MFillH fillh2(&mhn2, "", "FillH2");
419	MFillH fillh3(&mhn3, "", "FillH3");
420	MFillH fillh4(&mhn4, "", "FillH4");
421	MFillH filltp(&mhtp, "", "FillTriggerPos");
422	MFillH fillew(&mhew, "", "FillEvtWidth");
423	MFillH filled(&mhed, "", "FillMedDev");
424	fillh1.SetNameTab("Corsika", "Distribution as simulated by Corsika");
425	fillh2.SetNameTab("Detectable", "Distribution of events hit the detector");
426	fillh3.SetNameTab("Triggered", "Distribution of triggered events");
427	fillh4.SetNameTab("Cleaned", "Distribution after cleaning and cuts");
428	filltp.SetNameTab("TrigPos", "TriggerPosition w.r.t the first photon");
429	fillew.SetNameTab("EvtWidth", "Time between first and last photon hitting a detector");
430	filled.SetNameTab("MedDev", "Time between first and last photon hitting a detector");
431
432	MHPhotonEvent planeG(1, "HPhotonEventGround"); // Get from MaxImpact
433	MHPhotonEvent plane0(2, "HMirrorPlane0"); // Get from MReflector
434	//MHPhotonEvent plane1(2, "HMirrorPlane1");
435	MHPhotonEvent plane2(2, "HMirrorPlane2");
436	MHPhotonEvent plane3(2, "HMirrorPlane3");
437	MHPhotonEvent plane4(2, "HMirrorPlane4");
438	MHPhotonEvent planeF1(3, "HPhotonEventCamera"); // Get from MGeomCam
439	MHPhotonEvent planeF2(header.IsPointRun()?4:3, "HPhotonEventCones"); // Get from MGeomCam
440
441	plist.AddToList(&planeG);
442	plist.AddToList(&plane0);
443	plist.AddToList(&plane2);
444	plist.AddToList(&plane3);
445	plist.AddToList(&plane4);
446	plist.AddToList(&planeF1);
447	plist.AddToList(&planeF2);;
448
449	//MHPSF psf;
450
451	MFillH fillG(&planeG, "MPhotonEvent", "FillGround");
452	MFillH fill0(&plane0, "MirrorPlane0", "FillReflector");
453	//MFillH fill1(&plane1, "MirrorPlane1", "FillCamShadow");
454	MFillH fill2(&plane2, "MirrorPlane2", "FillCandidates");
455	MFillH fill3(&plane3, "MirrorPlane3", "FillReflected");
456	MFillH fill4(&plane4, "MirrorPlane4", "FillFocal");
457	MFillH fillF1(&planeF1, "MPhotonEvent", "FillCamera");
458	MFillH fillF2(&planeF2, "MPhotonEvent", "FillCones");
459	//MFillH fillP(&psf, "MPhotonEvent", "FillPSF");
460
461	fillG.SetNameTab("Ground", "Photon distribution at ground");
462	fill0.SetNameTab("Reflector", "Photon distribution at reflector plane w/o camera shadow");
463	//fill1.SetNameTab("CamShadow", "Photon distribution at reflector plane w/ camera shadow");
464	fill2.SetNameTab("Candidates", "'Can hit' photon distribution at reflector plane w/ camera shadow");
465	fill3.SetNameTab("Reflected", "Photon distribution after reflector projected to reflector plane");
466	fill4.SetNameTab("Focal", "Photon distribution at focal plane");
467	fillF1.SetNameTab("Camera", "Photon distribution which hit the detector");
468	fillF2.SetNameTab("Cones", "Photon distribution after cones");
469	//fillP.SetNameTab("PSF", "Photon distribution after cones for all mirrors");
470
471	// -------------------------------------------------------------------
472
473	// FIXME: Pedestal and Calibration runs should get P and C
474	const TString rule1(Form("s/cer([0-9]+)/%s\\/00$1_R_MonteCarlo.root/", Esc(fPathOut).Data()));
475	const TString rule2(Form("s/cer([0-9]+)/%s\\/00$1_Y_MonteCarlo.root/", Esc(fPathOut).Data()));
476	const TString rule3(Form("s/cer([0-9]+)/%s\\/00$1_D_MonteCarlo.root/", Esc(fPathOut).Data()));
477
478	MWriteRootFile write3a( 2, rule3, fOverwrite?"RECREATE":"NEW", "Camera file");
479	MWriteRootFile write3b( 2, rule3, fOverwrite?"RECREATE":"NEW", "Camera file");
480	MWriteRootFile write2a( 2, rule2, fOverwrite?"RECREATE":"NEW", "Signal file");
481	MWriteRootFile write2b( 2, rule2, fOverwrite?"RECREATE":"NEW", "Signal file");
482	MWriteRootFile write1a( 2, rule1, fOverwrite?"RECREATE":"NEW", "Reflector file");
483	MWriteRootFile write1b( 2, rule1, fOverwrite?"RECREATE":"NEW", "Reflector file");
484
485	write1a.SetName("WriteRefData");
486	write1b.SetName("WriteRefMC");
487	write2a.SetName("WriteSigData");
488	write2b.SetName("WriteSigMC");
489	write3a.SetName("WriteCamData");
490	write3b.SetName("WriteCamMC");
491
492	SetupCommonFileStructure(write1a);
493	SetupCommonFileStructure(write2a);
494	SetupCommonFileStructure(write3a);
495
496	// R: Dedicated file structureedicated events
497	write1a.AddContainer("MPhotonEvent", "Events");
498
499	// I: Dedicated file structureedicated events
500	write2a.AddContainer("MPedPhotFromExtractorRndm", "RunHeaders"); // FIXME: Needed for the signal files to be display in MARS
501	write2a.AddContainer("MSignalCam", "Events");
502
503	// D: Dedicated file structureedicated events
504	write3a.AddContainer("ElectronicNoise", "RunHeaders");
505	write3a.AddContainer("MRawEvtData", "Events");
506
507	// Basic MC data
508	write1b.AddContainer("MMcEvtBasic", "OriginalMC");
509	write2b.AddContainer("MMcEvtBasic", "OriginalMC");
510	write3b.AddContainer("MMcEvtBasic", "OriginalMC");
511
512	// -------------------------------------------------------------------
513
514	MGeomApply apply;
515
516	MSimPSF simpsf;
517
518	MSimGeomCam simgeom;
519	simgeom.SetNameGeomCam("GeomCones");
520
521	MSimCalibrationSignal simcal;
522	simcal.SetNameGeomCam("GeomCones");
523
524	// FIXME: Simulate photons before cones and QE!
525
526	MSimRandomPhotons simnsb;
527	simnsb.SetFreq(5.8, 0.004); // ph/m^2/nm/sr/ns NSB, 4MHz dark count rate
528	simnsb.SetNameGeomCam("GeomCones");
529
530	// FIXME: How do we handle star-light? We need the rate but we also
531	// need to process it through the mirror!
532
533	MSimAPD simapd;
534	simapd.SetNameGeomCam("GeomCones");
535
536	MSimExcessNoise simexcnoise;
537	MSimBundlePhotons simsum;
538	MSimSignalCam simsignal;
539	MSimCamera simcam;
540	MSimTrigger simtrig;
541	MSimReadout simdaq;
542
543	MContinue conttrig("TriggerPos.fVal<0", "ContTrigger");
544
545	MParameterD pulpos("IntendedPulsePos");
546	// FIXME: Set a default which could be 1/3 of the digitization window
547	// pulpos.SetVal(40); // [ns]
548	plist.AddToList(&pulpos);
549
550	MParameterD trigger("TriggerPos");
551	trigger.SetVal(0);
552	plist.AddToList(&trigger);
553
554	// -------------------------------------------------------------------
555
556	// Remove isolated pixels
557	MImgCleanStd clean(0, 0);
558	clean.SetCleanLvl0(0); // The level above which isolated pixels are kept
559	clean.SetCleanRings(0);
560	clean.SetMethod(MImgCleanStd::kAbsolute);
561
562	//clean.SetNamePedPhotCam("MPedPhotFromExtractorRndm");
563
564	MHillasCalc hcalc;
565	hcalc.Disable(MHillasCalc::kCalcConc);
566
567	MHCamEvent evt0a(/10/3, "Signal", "Average signal (absolute);;S [ph]");
568	MHCamEvent evt0c(/10/3, "MaxSignal", "Maximum signal (absolute);;S [ph]");
569	MHCamEvent evt0d(/11/8, "ArrTm", "Time after first photon;;T [ns]");
570	evt0a.SetErrorSpread(kFALSE);
571	evt0c.SetCollectMax();
572
573	MContinue cut("", "Cut");
574
575	MFillH fillx0a(&evt0a, "MSignalCam", "FillAvgSignal");
576	MFillH fillx0c(&evt0c, "MSignalCam", "FillMaxSignal");
577	MFillH fillx0d(&evt0d, "MSignalCam", "FillArrTm");
578	MFillH fillx1("MHHillas", "MHillas", "FillHillas");
579	MFillH fillx3("MHHillasSrc", "MHillasSrc", "FillHillasSrc");
580	MFillH fillx4("MHNewImagePar", "MNewImagePar", "FillNewImagePar");
581	MFillH fillth("MHThreshold", "", "FillThreshold");
582	MFillH fillca("MHCollectionArea", "", "FillTrigArea");
583
584	fillth.SetNameTab("Threshold");
585	fillca.SetNameTab("TrigArea");
586
587	// -------------------------------------------------------------------
588
589	if (header.IsDataRun())
590	{
591	tasks.AddToList(&read);
592	tasks.AddToList(&pointing);
593	tasks.AddToList(&simmmcs);
594	if (!fPathOut.IsNull() && !HasNullOut())
595	{
596	tasks.AddToList(&write1b);
597	tasks.AddToList(&write2b);
598	if (fCamera)
599	tasks.AddToList(&write3b);
600	}
601	// if (header.IsPointRun())
602	// tasks.AddToList(&stars);
603	if (1)
604	tasks.AddToList(&simatm); // Here because before fillh1
605	tasks.AddToList(&fillh1);
606	tasks.AddToList(&fillG);
607	if (!header.IsPointRun())
608	{
609	tasks.AddToList(&absapd);
610	tasks.AddToList(&absmir);
611	if (0)
612	tasks.AddToList(&simatm); // FASTER?
613	}
614	tasks.AddToList(&reflect);
615	if (!header.IsPointRun())
616	{
617	tasks.AddToList(&fill0);
618	//tasks.AddToList(&fill1);
619	tasks.AddToList(&fill2);
620	tasks.AddToList(&fill3);
621	tasks.AddToList(&fill4);
622	tasks.AddToList(&fillF1);
623	}
624	tasks.AddToList(&cones);
625	//if (header.IsPointRun())
626	// tasks.AddToList(&fillP);
627	tasks.AddToList(&cont1);
628	}
629	// -------------------------------
630	tasks.AddToList(&apply);
631	if (header.IsDataRun())
632	{
633	tasks.AddToList(&simpsf);
634	tasks.AddToList(&simgeom);
635	tasks.AddToList(&cont2);
636	if (!header.IsPointRun())
637	{
638	tasks.AddToList(&fillF2);
639	tasks.AddToList(&fillh2);
640	}
641	tasks.AddToList(&cont3);
642	}
643	if (fCamera)
644	{
645	if (header.IsPedestalRun() \|\| header.IsCalibrationRun())
646	tasks.AddToList(&simcal);
647	tasks.AddToList(&simnsb);
648	tasks.AddToList(&simapd);
649	tasks.AddToList(&simexcnoise);
650	}
651	tasks.AddToList(&simsum);
652	if (fCamera)
653	{
654	tasks.AddToList(&simcam);
655	if (header.IsDataRun() \|\| fForceTrigger)
656	tasks.AddToList(&simtrig);
657	tasks.AddToList(&conttrig);
658	tasks.AddToList(&simdaq);
659	}
660	tasks.AddToList(&simsignal); // What do we do if signal<0?
661	if (!fPathOut.IsNull() && !HasNullOut())
662	{
663	tasks.AddToList(&write1a);
664	if (!header.IsPedestalRun())
665	tasks.AddToList(&write2a);
666	if (fCamera)
667	tasks.AddToList(&write3a);
668	}
669	// -------------------------------
670	if (fCamera)
671	{
672	// FIXME: MHCollectionArea Trigger Area!
673	if (header.IsDataRun())
674	tasks.AddToList(&fillh3);
675	tasks.AddToList(&filltp);
676	}
677	if (header.IsDataRun())
678	{
679	tasks.AddToList(&fillew);
680	tasks.AddToList(&filled);
681	}
682	if (!header.IsPedestalRun())
683	{
684	tasks.AddToList(&fillx0a);
685	tasks.AddToList(&fillx0c);
686	if (!header.IsCalibrationRun())
687	tasks.AddToList(&clean);
688	tasks.AddToList(&hcalc);
689	tasks.AddToList(&cut);
690	tasks.AddToList(&fillx0d);
691	tasks.AddToList(&fillx1);
692	//tasks.AddToList(&fillx2);
693	tasks.AddToList(&fillx3);
694	tasks.AddToList(&fillx4);
695	//tasks.AddToList(&fillx5);
696	}
697	if (header.IsDataRun())
698	{
699	tasks.AddToList(&fillh4);
700	tasks.AddToList(&fillth);
701	tasks.AddToList(&fillca);
702	}
703	//-------------------------------------------
704
705	tasks.SetAccelerator(MTask::kAccDontReset\|MTask::kAccDontTime);
706
707	// Create and setup the eventloop
708	MEvtLoop evtloop(fName);
709	evtloop.SetParList(&plist);
710	evtloop.SetDisplay(fDisplay);
711	evtloop.SetLogStream(&gLog);
712	if (!SetupEnv(evtloop))
713	return kFALSE;
714
715	// FIXME: If pedestal file given use the values from this file
716	//-------------------------------------------
717
718	MGeomCam cam = static_cast<MGeomCam>(env2.GetCont());
719
720	if (binstr.IsDefault())
721	binstr.SetEdgesLin(150, -shape.GetWidth(),
722	header.GetFreqSampling()/1000.*header.GetNumSamples()+shape.GetWidth());
723
724	if (binsd.IsDefault() && cam)
725	binsd.SetEdgesLin(100, 0, cam->GetMaxRadius()*cam->GetConvMm2Deg());
726
727	if (binsd0.IsDefault() && cam)
728	binsd0.SetEdgesLin(100, 0, cam->GetMaxRadius()*cam->GetConvMm2Deg());
729
730
731	header.Print();
732
733	// FIXME: Display acceptance curves!
734
735	if (fDisplay)
736	{
737	TCanvas &c = fDisplay->AddTab("Info");
738	c.Divide(2,2);
739
740	c.cd(1);
741	gPad->SetBorderMode(0);
742	gPad->SetFrameBorderMode(0);
743	gPad->SetGridx();
744	gPad->SetGridy();
745	gROOT->SetSelectedPad(0);
746	shape.DrawClone()->SetBit(kCanDelete);
747
748	if (env0.GetCont() && (header.IsDataRun() \|\| header.IsPointRun()))
749	{
750	c.cd(3);
751	gPad->SetBorderMode(0);
752	gPad->SetFrameBorderMode(0);
753	gPad->SetGridx();
754	gPad->SetGridy();
755	gROOT->SetSelectedPad(0);
756	static_cast<MReflector*>(env0.GetCont())->DrawClone("line")->SetBit(kCanDelete);
757	}
758
759	if (fCamera)
760	{
761	if (env1.GetCont())
762	{
763	c.cd(2);
764	gPad->SetBorderMode(0);
765	gPad->SetFrameBorderMode(0);
766	gPad->SetGridx();
767	gPad->SetGridy();
768	gROOT->SetSelectedPad(0);
769	MHCamera c = new MHCamera(static_cast<MGeomCam&>(env1.GetCont()));
770	c->SetStats(kFALSE);
771	c->SetBit(MHCamera::kNoLegend);
772	c->SetBit(kCanDelete);
773	c->Draw();
774	}
775
776	if (cam)
777	{
778	c.cd(4);
779	gPad->SetBorderMode(0);
780	gPad->SetFrameBorderMode(0);
781	gPad->SetGridx();
782	gPad->SetGridy();
783	gROOT->SetSelectedPad(0);
784	MHCamera c = new MHCamera(cam);
785	c->SetStats(kFALSE);
786	c->SetBit(MHCamera::kNoLegend);
787	c->SetBit(kCanDelete);
788	c->Draw();
789	}
790	}
791
792	TCanvas &d = fDisplay->AddTab("Info2");
793	d.Divide(2,2);
794
795	d.cd(1);
796	gPad->SetBorderMode(0);
797	gPad->SetFrameBorderMode(0);
798	gPad->SetGridx();
799	gPad->SetGridy();
800	gROOT->SetSelectedPad(0);
801	splinepde.DrawClone()->SetBit(kCanDelete);
802
803	d.cd(2);
804	gPad->SetBorderMode(0);
805	gPad->SetFrameBorderMode(0);
806	gPad->SetGridx();
807	gPad->SetGridy();
808	gROOT->SetSelectedPad(0);
809	splinemirror.DrawClone()->SetBit(kCanDelete);
810
811	d.cd(3);
812	gPad->SetBorderMode(0);
813	gPad->SetFrameBorderMode(0);
814	gPad->SetGridx();
815	gPad->SetGridy();
816	gROOT->SetSelectedPad(0);
817	splinecones.DrawClone()->SetBit(kCanDelete);
818
819	d.cd(4);
820	gPad->SetBorderMode(0);
821	gPad->SetFrameBorderMode(0);
822	gPad->SetGridx();
823	gPad->SetGridy();
824	gROOT->SetSelectedPad(0);
825	MParSpline all = (MParSpline)splinepde.DrawClone();
826	//all->SetTitle("Combined acceptance");
827	all->SetBit(kCanDelete);
828	all->Multiply(*splinemirror.GetSpline());
829	}
830
831	//-------------------------------------------
832
833	// Execute first analysis
834	if (!evtloop.Eventloop(fMaxEvents))
835	{
836	gLog << err << GetDescriptor() << ": Failed." << endl;
837	return kFALSE;
838	}
839
840	//-------------------------------------------
841	// FIXME: Display Spline in tab
842	// FIXME: Display Reflector in tab
843	// FIXME: Display Routing(?) in tab
844	// FIXME: Display Camera(s) in tab
845	//-------------------------------------------
846
847	if (!WriteResult(plist))
848	return kFALSE;
849
850	*fLog << all << GetDescriptor() << ": Done." << endl << endl << endl;;
851
852	return kTRUE;
853	}

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

Download in other formats: