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

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