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

Last change on this file since 19713 was 19706, checked in by tbretz, 5 years ago
Improved naming convention.
File size: 49.3 KB
Line 
1/* ======================================================================== *\
2!
3! *
4! * This file is part of MARS, the MAGIC Analysis and Reconstruction
5! * Software. It is distributed to you in the hope that it can be a useful
6! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
7! * It is distributed WITHOUT ANY WARRANTY.
8! *
9! * Permission to use, copy, modify and distribute this software and its
10! * documentation for any purpose is hereby granted without fee,
11! * provided that the above copyright notice appear in all copies and
12! * that both that copyright notice and this permission notice appear
13! * in supporting documentation. It is provided "as is" without express
14! * or implied warranty.
15! *
16!
17!
18! Author(s): Thomas Bretz, 1/2009 <mailto:tbretz@astro.uni-wuerzburg.de>
19!
20! Copyright: MAGIC Software Development, 2000-2009
21!
22!
23\* ======================================================================== */
24
25/////////////////////////////////////////////////////////////////////////////
26//
27// MJSimulation
28//
29//
30// Force reading a corsika file even if the footer (RUNE-section) is missing
31// by setting fForceMode to kTRUE or from the resource file by
32//
33// ForceMode: Yes
34//
35//
36// To switch off the simulation of the camera electronics, use:
37//
38// Camera: Off
39//
40// Note, that the border of the camera and the propertied of the cones
41// are still simulated (simply because it is fast). Furthermore, this
42// switches off the trigger for the output, i.e. all data which deposits
43// at least one photon in at least one pixel is written to the output file.
44//
45//
46// In case of a pedestal or calibration run the artificial trigger can
47// be "switched off" and the cosmics trrigger "switched on" by setting
48// fForceTrigger to kTRUE or from the resource file by
49//
50// ForceTrigger: Yes
51//
52//
53/////////////////////////////////////////////////////////////////////////////
54#include "MJSimulation.h"
55
56#include <TEnv.h>
57#include <TCanvas.h>
58#include <iostream>
59
60// Core
61#include "MLog.h"
62#include "MLogManip.h"
63
64#include "MArgs.h"
65#include "MDirIter.h"
66#include "MParList.h"
67#include "MTaskList.h"
68#include "MEvtLoop.h"
69
70#include "MStatusDisplay.h"
71
72// Tasks
73#include "MCorsikaRead.h"
74#include "MContinue.h"
75#include "MFillH.h"
76#include "MGeomApply.h"
77#include "MParameterCalc.h"
78#include "MSrcPosCalc.h"
79#include "MHillasCalc.h"
80#include "MImgCleanStd.h"
81#include "MWriteRootFile.h"
82#include "MWriteFitsFile.h"
83
84#include "MSimMMCS.h"
85#include "MSimAbsorption.h"
86#include "MSimAtmosphere.h"
87#include "MSimReflector.h"
88#include "MSimPointingPos.h"
89#include "MSimPSF.h"
90#include "MSimGeomCam.h"
91#include "MSimSignalCam.h"
92#include "MSimAPD.h"
93#include "MSimExcessNoise.h"
94#include "MSimCamera.h"
95#include "MSimTrigger.h"
96#include "MSimReadout.h"
97#include "MSimRandomPhotons.h"
98#include "MSimBundlePhotons.h"
99#include "MSimCalibrationSignal.h"
100
101// Histograms
102#include "MBinning.h"
103
104#include "MHn.h"
105#include "MHCamera.h"
106#include "MHCamEvent.h"
107#include "MHPhotonEvent.h"
108
109// Container
110#include "MRawRunHeader.h"
111#include "MParameters.h"
112#include "MReflector.h"
113#include "MParEnv.h"
114#include "MSpline3.h"
115#include "MParSpline.h"
116#include "MGeomCam.h"
117#include "MPedestalCam.h"
118#include "MPedestalPix.h"
119
120ClassImp(MJSimulation);
121
122using namespace std;
123
124// --------------------------------------------------------------------------
125//
126// Default constructor.
127//
128// Sets fRuns to 0, fExtractor to NULL, fDataCheck to kFALSE
129//
130MJSimulation::MJSimulation(const char *name, const char *title)
131 : fForceMode(kFALSE), fCamera(kTRUE), fForceTrigger(kFALSE),
132 fWriteFitsFile(kFALSE), fOperationMode(kModeData), fRunNumber(-1)
133{
134 fName = name ? name : "MJSimulation";
135 fTitle = title ? title : "Standard analysis and reconstruction";
136}
137
138Bool_t MJSimulation::CheckEnvLocal()
139{
140 fForceMode = GetEnv("ForceMode", fForceMode);
141 fForceTrigger = GetEnv("ForceTrigger", fForceTrigger);
142 fCamera = GetEnv("Camera", fCamera);
143
144 return kTRUE;
145}
146
147Bool_t MJSimulation::WriteResult(const MParList &plist, const MSequence &seq, Int_t num)
148{
149 if (fPathOut.IsNull())
150 {
151 *fLog << inf << "No output path specified via SetPathOut - no output written." << endl;
152 return kTRUE;
153 }
154
155 TObjArray cont;
156 cont.Add(const_cast<TEnv*>(GetEnv()));
157 if (seq.IsValid())
158 cont.Add(const_cast<MSequence*>(&seq));
159
160 cont.Add(plist.FindObject("PulseShape"));
161 cont.Add(plist.FindObject("Reflector"));
162 cont.Add(plist.FindObject("MGeomCam"));
163 cont.Add(plist.FindObject("GeomCones"));
164
165 TNamed cmdline("CommandLine", fCommandLine.Data());
166 cont.Add(&cmdline);
167
168 if (fDisplay)
169 {
170 TString title = "-- Ceres";
171 if (seq.IsValid())
172 {
173 title += ": ";
174 title += seq.GetSequence();
175 }
176 title += " --";
177 fDisplay->SetTitle("Ceres", kFALSE);
178
179 cont.Add(fDisplay);
180 }
181
182 const TString name = Form("ceres%08d.root", num);
183 return WriteContainer(cont, name, "RECREATE");
184}
185
186void MJSimulation::SetupHist(MHn &hist) const
187{
188 hist.AddHist("MCorsikaEvtHeader.fTotalEnergy");
189 hist.InitName("Energy");
190 hist.InitTitle("Energy;E [GeV]");
191 hist.SetLog(kTRUE, kTRUE, kFALSE);
192
193 hist.AddHist("MPhotonEvent.GetNumExternal");
194 hist.InitName("Size");
195 hist.InitTitle("Size;S [#]");
196 hist.SetLog(kTRUE, kTRUE, kFALSE);
197
198 /*
199 hist.AddHist("-MCorsikaEvtHeader.fX/100","-MCorsikaEvtHeader.fY/100");
200 hist.SetDrawOption("colz");
201 hist.InitName("Impact;Impact;Impact");
202 hist.InitTitle("Impact;West <--> East [m];South <--> North [m]");
203 hist.SetAutoRange();
204 */
205
206 hist.AddHist("MCorsikaEvtHeader.GetImpact/100");
207 hist.InitName("Impact");
208 hist.InitTitle("Impact;Impact [m]");
209 hist.SetAutoRange();
210
211 hist.AddHist("MCorsikaEvtHeader.fFirstInteractionHeight/100000");
212 hist.InitName("Height");
213 hist.InitTitle("FirstInteractionHeight;h [km]");
214
215 hist.AddHist("(MCorsikaEvtHeader.fAz+MCorsikaRunHeader.fMagneticFieldAz)*TMath::RadToDeg()", "MCorsikaEvtHeader.fZd*TMath::RadToDeg()");
216 hist.InitName("SkyOrigin;Az;Zd");
217 hist.InitTitle("Sky Origin;Az [#circ];Zd [#circ]");
218 hist.SetDrawOption("colz");
219 hist.SetAutoRange();
220
221 hist.AddHist("IncidentAngle.fVal");
222 hist.InitName("ViewCone");
223 hist.InitTitle("Incident Angle;#alpha [#circ]");
224}
225
226template<class T>
227void MJSimulation::SetupCommonFileStructure(T &write) const
228{
229 // Common run headers
230 write.AddContainer("MMcCorsikaRunHeader", "RunHeaders", kFALSE, 1);
231 write.AddContainer("MCorsikaRunHeader", "RunHeaders", kFALSE, 1);
232 write.AddContainer("MRawRunHeader", "RunHeaders", kTRUE, 1);
233 write.AddContainer("MGeomCam", "RunHeaders", kTRUE, 1);
234 write.AddContainer("MMcRunHeader", "RunHeaders", kTRUE, 1);
235
236 // Common events
237 write.AddContainer("MCorsikaEvtHeader", "Events", kFALSE);
238 write.AddContainer("MRawEvtHeader", "Events");
239 write.AddContainer("MMcEvt", "Events", kFALSE);
240 write.AddContainer("MMcEvtBasic", "Events", kFALSE);
241 write.AddContainer("IncidentAngle", "Events", kFALSE);
242 write.AddContainer("MPointingPos", "Events", kFALSE);
243 write.AddContainer("MSimSourcePos", "Events", kFALSE);
244 write.AddContainer("MSrcPosCam", "Events", kFALSE);
245}
246
247void 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
323void 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
335void 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
363void 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","BinningWidth");
388 plist.FindCreateObj("MBinning","BinningLength");
389 plist.FindCreateObj("MBinning","BinningDistC");
390 plist.FindCreateObj("MBinning","BinningDist");
391 plist.FindCreateObj("MBinning","BinningTrigPos");
392}
393
394Bool_t MJSimulation::Process(const MArgs &args, const MSequence &seq)
395{
396 // --------------------------------------------------------------------------------
397 // --------------------------------------------------------------------------------
398 // Initialization
399 // --------------------------------------------------------------------------------
400 // --------------------------------------------------------------------------------
401
402 // --------------------------------------------------------------------------------
403 // Logging output:
404 // --------------------------------------------------------------------------------
405 // - description of the job itself
406 // - list of the to processed sequence
407 *fLog << inf;
408 fLog->Separator(GetDescriptor());
409
410 if (!CheckEnv())
411 return kFALSE;
412
413 if (seq.IsValid())
414 *fLog << fSequence.GetFileName() << endl;
415 else
416 *fLog << "Input: " << args.GetNumArguments() << "-files" << endl;
417 *fLog << endl;
418
419 MDirIter iter;
420 if (seq.IsValid() && seq.GetRuns(iter, MSequence::kCorsika)<=0)
421 {
422 *fLog << err << "ERROR - Sequence valid but without files." << endl;
423 return kFALSE;
424 }
425 // --------------------------------------------------------------------------------
426 // Setup MParList and MTasklist
427 // --------------------------------------------------------------------------------
428 MParList plist;
429 plist.AddToList(this); // take care of fDisplay!
430 // setup TaskList
431 MTaskList tasks;
432 plist.AddToList(&tasks);
433
434 // --------------------------------------------------------------------------------
435 // --------------------------------------------------------------------------------
436 // Parameter Container Setup
437 // --------------------------------------------------------------------------------
438 // --------------------------------------------------------------------------------
439
440 // --------------------------------------------------------------------------------
441 // Setup container for the reflector, the cones and the camera geometry
442 // --------------------------------------------------------------------------------
443 // FIXME: Allow empty GeomCones!
444 MParEnv env0("Reflector");
445 MParEnv env1("GeomCones");
446 MParEnv env2("MGeomCam");
447 env0.SetClassName("MOptics");
448 env1.SetClassName("MGeomCam");
449 env2.SetClassName("MGeomCam");
450 plist.AddToList(&env0);
451 plist.AddToList(&env1);
452 plist.AddToList(&env2);
453 // --------------------------------------------------------------------------------
454 // Setup container for the "camera array" of the extracted number of photons
455 // --------------------------------------------------------------------------------
456 // from ExtractorRndm
457 plist.FindCreateObj("MPedPhotCam", "MPedPhotFromExtractorRndm");
458 // --------------------------------------------------------------------------------
459 // Setup spline containers for:
460 // --------------------------------------------------------------------------------
461 // - the pulse shape
462 // - the different photon acceptance splines
463 MParSpline shape("PulseShape");
464 MParSpline splinepde("PhotonDetectionEfficiency");
465 MParSpline splinemirror("MirrorReflectivity");
466 MParSpline splinecones("ConesAngularAcceptance");
467 MParSpline splinecones2("ConesTransmission");
468 MParSpline splineAdditionalPhotonAcceptance("AdditionalPhotonAcceptance");
469 plist.AddToList(&shape);
470 plist.AddToList(&splinepde);
471 plist.AddToList(&splinemirror);
472 plist.AddToList(&splinecones);
473 plist.AddToList(&splinecones2);
474 plist.AddToList(&splineAdditionalPhotonAcceptance);
475
476 // --------------------------------------------------------------------------------
477 // Setup container for the MC run header and the Raw run header
478 // --------------------------------------------------------------------------------
479 plist.FindCreateObj("MMcRunHeader");
480
481 MRawRunHeader header;
482 header.SetValidMagicNumber();
483 SetupHeaderOperationMode(header);
484 // FIXME: Move to MSimPointingPos, MSimCalibrationSignal
485 // Can we use this as input for MSimPointingPos?
486 header.SetObservation("On", "MonteCarlo");
487 plist.AddToList(&header);
488
489 // --------------------------------------------------------------------------------
490 // Setup container for the intended pulse position and for the trigger position
491 // --------------------------------------------------------------------------------
492 MParameterD pulpos("IntendedPulsePos");
493 // FIXME: Set a default which could be 1/3 of the digitization window
494 // pulpos.SetVal(40); // [ns]
495 plist.AddToList(&pulpos);
496
497 MParameterD trigger("TriggerPos");
498 trigger.SetVal(0);
499 plist.AddToList(&trigger);
500
501
502 // --------------------------------------------------------------------------------
503 // Setup container for the fix time offset, for residual time spread and for gapd time jitter
504 // --------------------------------------------------------------------------------
505 // Dominik and Sebastian on: fix time offsets
506 MMatrix fix_time_offsets_between_pixels_in_ns("MFixTimeOffset", 0);
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 if (read.GetNumFiles()==0 && fOperationMode==kModeData)
544 {
545 *fLog << err << "No files to read." << endl;
546 return kFALSE;
547 }
548
549 // --------------------------------------------------------------------------------
550 // Precut
551 // --------------------------------------------------------------------------------
552 MContinue precut("", "PreCut");
553 precut.IsInverted();
554 precut.SetAllowEmpty();
555
556 // --------------------------------------------------------------------------------
557 // MSimMMCS (don't know what this is doing) and calculation of the incident angle
558 // --------------------------------------------------------------------------------
559 MSimMMCS simmmcs;
560 const TString sin2 = "sin(MCorsikaEvtHeader.fZd)*sin(MCorsikaRunHeader.fZdMin*TMath::DegToRad())";
561 const TString cos2 = "cos(MCorsikaEvtHeader.fZd)*cos(MCorsikaRunHeader.fZdMin*TMath::DegToRad())";
562 const TString cos = "cos(MCorsikaEvtHeader.fAz-MCorsikaRunHeader.fAzMin*TMath::DegToRad())";
563
564 const TString form = "acos("+sin2+"*"+cos+"+"+cos2+")*TMath::RadToDeg()";
565
566 MParameterCalc calcangle(form, "CalcIncidentAngle");
567 calcangle.SetNameParameter("IncidentAngle");
568
569 // --------------------------------------------------------------------------------
570 // Absorption tasks
571 // --------------------------------------------------------------------------------
572 // - Atmosphere (simatm)
573 // - photon detection efficiency (absapd)
574 // - mirror reflectivity (absmir)
575 // - angular acceptance of winston cones (cones)
576 // - transmission of winston cones (cones2)
577 // - additional photon acceptance (only motivated, not measured) (additionalPhotonAcceptance)
578 MSimAtmosphere simatm;
579 MSimAbsorption absapd("SimPhotonDetectionEfficiency");
580 MSimAbsorption absmir("SimMirrorReflectivity");
581 MSimAbsorption cones("SimConesAngularAcceptance");
582 MSimAbsorption cones2("SimConesTransmission");
583 MSimAbsorption additionalPhotonAcceptance("SimAdditionalPhotonAcceptance");
584 absapd.SetParName("PhotonDetectionEfficiency");
585 absmir.SetParName("MirrorReflectivity");
586 cones.SetParName("ConesAngularAcceptance");
587 cones.SetUseTheta();
588 cones2.SetParName("ConesTransmission");
589 additionalPhotonAcceptance.SetParName("AdditionalPhotonAcceptance");
590 additionalPhotonAcceptance.SetForce(kTRUE);
591
592 // --------------------------------------------------------------------------------
593 // Simulating pointing position and simulating reflector
594 // --------------------------------------------------------------------------------
595 MSimPointingPos pointing;
596 MSrcPosCalc srcposcam;
597
598 MSimReflector reflect;
599 reflect.SetNameGeomCam("GeomCones");
600 reflect.SetNameReflector("Reflector");
601// MSimStarField stars;
602
603 // --------------------------------------------------------------------------------
604 // Continue tasks
605 // --------------------------------------------------------------------------------
606 // - Processing of the current events stops, if there aren't at least two photons (see cont3)
607 MContinue cont1("MPhotonEvent.GetNumPhotons<1", "ContEmpty1");
608 MContinue cont2("MPhotonEvent.GetNumPhotons<1", "ContEmpty2");
609 MContinue cont3("MPhotonEvent.GetNumPhotons<2", "ContEmpty3");
610 cont1.SetAllowEmpty();
611 cont2.SetAllowEmpty();
612 cont3.SetAllowEmpty();
613
614 // --------------------------------------------------------------------------------
615 // --------------------------------------------------------------------------------
616 // Tasks Setup (Cones, SiPM-Simulation, RandomPhotons, Camera, Trigger, DAQ)
617 // --------------------------------------------------------------------------------
618 // --------------------------------------------------------------------------------
619
620 // --------------------------------------------------------------------------------
621 // GeomApply, SimPSF, SimGeomCam, SimCalibrationSignal
622 // --------------------------------------------------------------------------------
623 // - MGeomApply only resizes all MGeomCam parameter container to the actual
624 // number of pixels
625 // - MSimPSF applies an additional smearing of the photons on the camera window
626 // on default it is switched of (set to 0), but can be applied by setting:
627 // MSimPSF.fSigma: <value>
628 // in the config file.
629 // Be aware, that there is also a smearing of the photons implemented in
630 // MSimReflector by setting:
631 // MReflector.SetSigmaPSF: <value>
632 // in the config file. This is at the moment done in the default config file.
633 // - MSimGeomCam identifies which pixel was hit by which photon
634 // - MSimCalibrationSignal is only used for simulated calibration runs
635 MGeomApply apply;
636
637 MSimPSF simpsf;
638
639 MSimGeomCam simgeom;
640 simgeom.SetNameGeomCam("GeomCones");
641
642 MSimCalibrationSignal simcal;
643 simcal.SetNameGeomCam("GeomCones");
644
645 // --------------------------------------------------------------------------------
646 // Simulation of random photons
647 // --------------------------------------------------------------------------------
648 // - be aware that here default values for the nsb rate and for the dark
649 // count rate are set. They will be overwritten if the values are defined
650 // in the config file
651 // - if you want to simulate a specific nsb rate you have to set the filename
652 // in the config file to an empty string and then you can specify the NSB rate:
653 // MSimRandomPhotons.FileNameNSB:
654 // MSimRandomPhotons.FrequencyNSB: 0.0
655
656 // Sky Quality Meter: 21.82M = 2.02e-4cd/m^2
657 // 1cd = 12.566 lm / 4pi sr
658 // FIXME: Simulate photons before cones and QE!
659 MSimRandomPhotons simnsb;
660 simnsb.SetFreq(5.8, 0.004); // ph/m^2/nm/sr/ns NSB, 4MHz dark count rate
661 simnsb.SetNameGeomCam("GeomCones");
662 // FIXME: How do we handle star-light? We need the rate but we also
663 // need to process it through the mirror!
664
665 // --------------------------------------------------------------------------------
666 // Simulation from the SiPM to the DAQ
667 // --------------------------------------------------------------------------------
668 // - MSimAPD simulates the whole behaviour of the SiPMs:
669 // (dead time of cells, crosstalk, afterpulses)
670 // - MSimExcessNoise adds a spread on the weight of each signal in the SiPMs
671 // - MSimBundlePhotons restructured the photon list of MPhotonEvent via a look
672 // up table. At the moment the tasks is skipped, due to the fact that there
673 // is no file name specified in the config file
674 // - MSimSignalCam only fills a SignalCam container for the cherenkov photons
675 // - MSimCamera simulates the behaviour of the camera:
676 // (electronic noise, accoupling, time smearing and offset for the photons,
677 // pulses injected by all photons)
678 // - MSimTrigger simulates the behaviour of the Trigger:
679 // (Adding patch voltages, applying clipping, applying discriminator, applying
680 // time over threshold and a possible trigger logic)
681 // - MContinue conttrig stops the processing of the current event if the trigger
682 // position is negativ (and therefore not valid)
683 // - MSimReadout simulates the behaviour of the readout:
684 // (Digitization and saturation of the ADC)
685 MSimAPD simapd;
686 simapd.SetNameGeomCam("GeomCones");
687
688 MSimExcessNoise simexcnoise;
689 MSimBundlePhotons simsum;
690 MSimSignalCam simsignal;
691 MSimCamera simcam;
692 MSimTrigger simtrig;
693 MContinue conttrig("TriggerPos.fVal<0", "ContTrigger");
694 MSimReadout simdaq;
695
696 // --------------------------------------------------------------------------------
697 // Standard analysis with hillas parameters for the true cherenkov photons
698 // --------------------------------------------------------------------------------
699 // Remove isolated pixels
700 MImgCleanStd clean(0, 0);
701 //clean.SetCleanLvl0(0); // The level above which isolated pixels are kept
702 clean.SetCleanRings(0);
703 clean.SetMethod(MImgCleanStd::kAbsolute);
704
705 //clean.SetNamePedPhotCam("MPedPhotFromExtractorRndm");
706
707 MHillasCalc hcalc;
708 hcalc.Disable(MHillasCalc::kCalcConc);
709
710
711 // --------------------------------------------------------------------------------
712 // Setup of histogram containers and the corresponding fill tasks
713 // --------------------------------------------------------------------------------
714 // Here is no functionality for the simulation, it is only visualization via
715 // showplot
716 MHn mhn1, mhn2, mhn3, mhn4;
717 SetupHist(mhn1);
718 SetupHist(mhn2);
719 SetupHist(mhn3);
720 SetupHist(mhn4);
721
722 MH3 mhtp("TriggerPos.fVal-IntendedPulsePos.fVal-PulseShape.GetWidth*1000/MRawRunHeader.GetFreqSampling");
723 mhtp.SetName("TrigPos");
724 mhtp.SetTitle("Trigger position w.r.t. the first photon (=Pulse.Xmin) hitting a detector [ns];#Delta T [ns]");
725
726 MH3 mhew("MPhotonStatistics.fLength");
727 mhew.SetName("TotLength");
728 mhew.SetTitle("Time between first and last photon hitting a detector;#Delta T [ns]");
729
730 MH3 mhed("MPhotonStatistics.fTimeMedDev");
731 mhed.SetName("MedLength");
732 mhed.SetTitle("Median deviation (1#sigma);#sigma(#Delta T) [ns]");
733
734 MFillH fillh1(&mhn1, "", "FillCorsika");
735 MFillH fillh2(&mhn2, "", "FillH2");
736 MFillH fillh3(&mhn3, "", "FillH3");
737 MFillH fillh4(&mhn4, "", "FillH4");
738 MFillH filltp(&mhtp, "", "FillTriggerPos");
739 MFillH fillew(&mhew, "", "FillEvtWidth");
740 MFillH filled(&mhed, "", "FillMedDev");
741 fillh1.SetNameTab("Corsika", "Distribution as simulated by Corsika");
742 fillh2.SetNameTab("Detectable", "Distribution of events hit the detector");
743 fillh3.SetNameTab("Triggered", "Distribution of triggered events");
744 fillh4.SetNameTab("Cleaned", "Distribution after cleaning and cuts");
745 filltp.SetNameTab("TrigPos", "TriggerPosition w.r.t the first photon");
746 fillew.SetNameTab("EvtWidth", "Time between first and last photon hitting a detector");
747 filled.SetNameTab("MedDev", "Median deviation of arrival time of photons hitting a detector");
748
749 MHPhotonEvent planeG(1, "HPhotonEventGround"); // Get from MaxImpact
750 MHPhotonEvent plane0(2, "HMirrorPlane0"); // Get from MReflector
751 //MHPhotonEvent plane1(2, "HMirrorPlane1");
752 MHPhotonEvent plane2(2, "HMirrorPlane2");
753 MHPhotonEvent plane3(2, "HMirrorPlane3");
754 MHPhotonEvent plane4(2, "HMirrorPlane4");
755 MHPhotonEvent planeF1(6, "HPhotonEventCamera"); // Get from MGeomCam
756 MHPhotonEvent planeF2(header.IsPointRun()?4:6, "HPhotonEventCones"); // Get from MGeomCam
757
758 plist.AddToList(&planeG);
759 plist.AddToList(&plane0);
760 plist.AddToList(&plane2);
761 plist.AddToList(&plane3);
762 plist.AddToList(&plane4);
763 plist.AddToList(&planeF1);
764 plist.AddToList(&planeF2);
765
766 //MHPSF psf;
767
768 MFillH fillG(&planeG, "MPhotonEvent", "FillGround");
769 MFillH fill0(&plane0, "MirrorPlane0", "FillReflector");
770 //MFillH fill1(&plane1, "MirrorPlane1", "FillCamShadow");
771 MFillH fill2(&plane2, "MirrorPlane2", "FillCandidates");
772 MFillH fill3(&plane3, "MirrorPlane3", "FillReflected");
773 MFillH fill4(&plane4, "MirrorPlane4", "FillFocal");
774 MFillH fillF1(&planeF1, "MPhotonEvent", "FillCamera");
775 MFillH fillF2(&planeF2, "MPhotonEvent", "FillCones");
776 //MFillH fillP(&psf, "MPhotonEvent", "FillPSF");
777
778 fillG.SetNameTab("Ground", "Photon distribution at ground");
779 fill0.SetNameTab("Reflector", "Photon distribution at reflector plane w/o camera shadow");
780 //fill1.SetNameTab("CamShadow", "Photon distribution at reflector plane w/ camera shadow");
781 fill2.SetNameTab("Candidates", "*Can hit* photon distribution at reflector plane w/ camera shadow");
782 fill3.SetNameTab("Reflected", "Photon distribution after reflector projected to reflector plane");
783 fill4.SetNameTab("Focal", "Photon distribution at focal plane");
784 fillF1.SetNameTab("Camera", "Photon distribution which hit the detector");
785 fillF2.SetNameTab("Cones", "Photon distribution after cones");
786 //fillP.SetNameTab("PSF", "Photon distribution after cones for all mirrors");
787
788
789 MHCamEvent evt0a(/*10*/3, "Signal", "Average signal (absolute);;S [ph]");
790 MHCamEvent evt0c(/*10*/3, "MaxSignal", "Maximum signal (absolute);;S [ph]");
791 MHCamEvent evt0d(/*11*/8, "ArrTm", "Time after first photon;;T [ns]");
792 evt0a.SetErrorSpread(kFALSE);
793 evt0c.SetCollectMax();
794
795 MContinue cut("", "Cut");
796
797 MFillH fillx0a(&evt0a, "MSignalCam", "FillAvgSignal");
798 MFillH fillx0c(&evt0c, "MSignalCam", "FillMaxSignal");
799 MFillH fillx0d(&evt0d, "MSignalCam", "FillArrTm");
800 MFillH fillx1("MHHillas", "MHillas", "FillHillas");
801 MFillH fillx3("MHHillasSrc", "MHillasSrc", "FillHillasSrc");
802 MFillH fillx4("MHNewImagePar", "MNewImagePar", "FillNewImagePar");
803 MFillH fillth("MHThreshold", "", "FillThreshold");
804 MFillH fillca("MHCollectionArea", "", "FillTrigArea");
805
806 fillth.SetNameTab("Threshold");
807 fillca.SetNameTab("TrigArea");
808
809 // --------------------------------------------------------------------------------
810 // Formating of the outputfilepath
811 // --------------------------------------------------------------------------------
812 if (!fFileOut.IsNull())
813 {
814 const Ssiz_t dot = fFileOut.Last('.');
815 const Ssiz_t slash = fFileOut.Last('/');
816 if (dot>slash)
817 fFileOut = fFileOut.Remove(dot);
818 }
819
820 const char *regex = fRunNumber<0 ?
821 "s/[cC][eE][rR]([0-9]+)(-[tT][eE][lL]([0-9]+))?/%s\\/00$1.%03d%%s_MonteCarlo$3.root/" :
822 "s/[cC][eE][rR][0-9]+([0-9][0-9][0-9])(-[tT][eE][lL]([0-9]+))?/%s\\/%08d.$1%%s_MonteCarlo$3.root/";
823
824 const char *fmt = fFileOut.IsNull() ?
825 Form(regex, Esc(fPathOut).Data(), header.GetRunNumber()) :
826 Form("%s/%s%%s.root", Esc(fPathOut).Data(), Esc(fFileOut).Data());
827
828 const TString rule1(Form(fmt, fFileOut.IsNull()?"_R":""));
829 const TString rule2(Form(fmt, "_Y"));
830 const TString rule4(Form(fmt, "_I"));
831 TString rule3(Form(fmt, Form("_%c", header.GetRunTypeChar())));
832
833 // --------------------------------------------------------------------------------
834 // Setup of the outputfile tasks
835 // --------------------------------------------------------------------------------
836 MWriteRootFile write4a( 2, rule4, fOverwrite?"RECREATE":"NEW", "Star file");
837 MWriteRootFile write4b( 2, rule4, fOverwrite?"RECREATE":"NEW", "Star file");
838 MWriteRootFile write3b( 2, rule3, fOverwrite?"RECREATE":"NEW", "Camera file");
839 MWriteRootFile write2a( 2, rule2, fOverwrite?"RECREATE":"NEW", "Signal file");
840 MWriteRootFile write2b( 2, rule2, fOverwrite?"RECREATE":"NEW", "Signal file");
841 MWriteRootFile write1a( 2, rule1, fOverwrite?"RECREATE":"NEW", "Reflector file");
842 MWriteRootFile write1b( 2, rule1, fOverwrite?"RECREATE":"NEW", "Reflector file");
843
844 if (fWriteFitsFile)
845 rule3.ReplaceAll(".root", ".fits");
846
847 MWriteFitsFile write3af( 2, rule3, fOverwrite?"RECREATE":"NEW", "Camera file");
848
849 if (fWriteFitsFile)
850 rule3 = "/dev/null";
851
852 MWriteRootFile write3ar( 2, rule3, fOverwrite?"RECREATE":"NEW", "Camera file");
853
854 MTask &write3a = fWriteFitsFile ? static_cast<MTask&>(write3af) : static_cast<MTask&>(write3ar);
855
856 //SetupHeaderKeys(write3af,header);
857 SetupVetoColumns(write3af);
858
859 write3af.SetBytesPerSample("Data", 2);
860
861 write1a.SetName("WriteRefData");
862 write1b.SetName("WriteRefMC");
863 write2a.SetName("WriteSigData");
864 write2b.SetName("WriteSigMC");
865 write3a.SetName("WriteCamData");
866 write3b.SetName("WriteCamMC");
867 write4a.SetName("WriteImgData");
868 write4b.SetName("WriteImgMC");
869
870 SetupCommonFileStructure(write1a);
871 SetupCommonFileStructure(write2a);
872 SetupCommonFileStructure(write3ar);
873 SetupCommonFileStructure(write3af);
874 SetupCommonFileStructure(write4a);
875
876 // R: Dedicated file structure
877 write1a.AddContainer("MPhotonEvent", "Events");
878
879 // Y: Dedicated file structure
880 write2a.AddContainer("MPedPhotFromExtractorRndm", "RunHeaders", kTRUE, 1); // FIXME: Needed for the signal files to be display in MARS
881 write2a.AddContainer("MSignalCam", "Events");
882
883 // D: Dedicated file structure
884 write3af.AddContainer("ElectronicNoise", "RunHeaders", kTRUE, 1);
885 write3af.AddContainer("IntendedPulsePos", "RunHeaders", kTRUE, 1);
886 write3af.AddContainer("ResidualTimeSpread", "RunHeaders", kTRUE, 1);
887 write3af.AddContainer("GapdTimeJitter", "RunHeaders", kTRUE, 1);
888 write3af.AddContainer("MRawEvtData", "Events");
889 write3af.AddContainer("MTruePhotonsPerPixelCont", "Events");
890
891 write3ar.AddContainer("ElectronicNoise", "RunHeaders", kTRUE, 1);
892 write3ar.AddContainer("IntendedPulsePos", "RunHeaders", kTRUE, 1);
893 write3ar.AddContainer("MRawEvtData", "Events");
894 // It doesn't make much sene to write this information
895 // to the file because the units are internal units not
896 // related to the output samples
897 // if (header.IsDataRun() || fForceTrigger)
898 // write3a.AddContainer("TriggerPos", "Events");
899
900 // I: Dedicated file structure
901 write4a.AddContainer("MHillas", "Events");
902 write4a.AddContainer("MHillasSrc", "Events");
903 write4a.AddContainer("MImagePar", "Events");
904 write4a.AddContainer("MNewImagePar", "Events");
905
906 // Basic MC data
907 write1b.AddContainer("MMcEvtBasic", "OriginalMC");
908 write2b.AddContainer("MMcEvtBasic", "OriginalMC");
909 write3b.AddContainer("MMcEvtBasic", "OriginalMC");
910 write4b.AddContainer("MMcEvtBasic", "OriginalMC");
911
912
913 // --------------------------------------------------------------------------------
914 // --------------------------------------------------------------------------------
915 // Filling of tasks in tasklist
916 // --------------------------------------------------------------------------------
917 // --------------------------------------------------------------------------------
918
919 if (header.IsDataRun())
920 {
921 tasks.AddToList(&read); // Reading Corsika
922 tasks.AddToList(&precut); // Precut
923 tasks.AddToList(&pointing); // Simulating pointing
924 tasks.AddToList(&srcposcam); // calculate origin in camera
925 tasks.AddToList(&simmmcs); // Simulating MMCS
926 if (!fPathOut.IsNull() && !HasNullOut()) // Write Tasks for corsika infos
927 {
928 //tasks.AddToList(&write1b);
929 tasks.AddToList(&write2b);
930 if (fCamera)
931 tasks.AddToList(&write3b);
932 if (header.IsDataRun())
933 tasks.AddToList(&write4b);
934 }
935 // if (header.IsPointRun())
936 // tasks.AddToList(&stars);
937 if (1)
938 tasks.AddToList(&simatm); // Here because before fillh1 // atmosphere absorption
939 tasks.AddToList(&calcangle); // calculation incident angle
940 tasks.AddToList(&fillh1); // fill histogram task
941 tasks.AddToList(&fillG); // fill histogram task
942 if (!header.IsPointRun())
943 {
944 tasks.AddToList(&absapd); // photon detection efficiency of the apds
945 tasks.AddToList(&absmir); // mirror reflectivity
946 tasks.AddToList(&additionalPhotonAcceptance); // addition photon acceptance
947 if (0)
948 tasks.AddToList(&simatm); // FASTER? // be aware this is 'commented'
949 }
950 tasks.AddToList(&reflect); // Simulation of the reflector
951 if (!header.IsPointRun())
952 {
953 tasks.AddToList(&fill0); // fill histogram task
954 //tasks.AddToList(&fill1);
955 tasks.AddToList(&fill2); // fill histogram task
956 tasks.AddToList(&fill3); // fill histogram task
957 tasks.AddToList(&fill4); // fill histogram task
958 tasks.AddToList(&fillF1); // fill histogram task
959 }
960 tasks.AddToList(&cones); // angular acceptance of winston cones
961 tasks.AddToList(&cones2); // transmission of winston cones
962 //if (header.IsPointRun())
963 // tasks.AddToList(&fillP);
964 tasks.AddToList(&cont1); // continue if at least 2 photons
965 }
966 // -------------------------------
967 tasks.AddToList(&apply); // apply geometry to MGeomCam containers
968 if (header.IsDataRun())
969 {
970 tasks.AddToList(&simpsf); // simulates additional psf (NOT used in default mode)
971 tasks.AddToList(&simgeom); // tag which pixel is hit by which photon
972 tasks.AddToList(&cont2); // continue if at least 2 photons
973 if (!header.IsPointRun())
974 {
975 tasks.AddToList(&fillF2); // fill histogram task
976 tasks.AddToList(&fillh2); // fill histogram task
977 }
978 tasks.AddToList(&cont3); // continue if at least 3 photons
979 }
980 if (fCamera)
981 {
982 if (header.IsPedestalRun() || header.IsCalibrationRun())
983 tasks.AddToList(&simcal); // add calibration signal for calibration runs
984 tasks.AddToList(&simnsb); // simulate nsb and dark counts
985 tasks.AddToList(&simapd); // simulate SiPM behaviour (dead time, crosstalk ...)
986 tasks.AddToList(&simexcnoise); // add excess noise
987 }
988 tasks.AddToList(&simsum); // bundle photons (NOT used in default mode)
989 if (fCamera)
990 {
991 tasks.AddToList(&simcam); // simulate camera behaviour (creates analog signal)
992 if (header.IsDataRun() || fForceTrigger)
993 tasks.AddToList(&simtrig); // simulate trigger
994 tasks.AddToList(&conttrig); // continue if trigger pos is valid
995 tasks.AddToList(&simdaq); // simulate data aquisition
996 }
997 tasks.AddToList(&simsignal); // What do we do if signal<0? // fill MSimSignal container
998 if (!fPathOut.IsNull() && !HasNullOut())
999 {
1000 //tasks.AddToList(&write1a);
1001 if (!header.IsPedestalRun())
1002 tasks.AddToList(&write2a); // outputtask
1003 if (fCamera)
1004 tasks.AddToList(&write3a); // outputtask (this is the raw data writing tasks)
1005 }
1006 // -------------------------------
1007 if (fCamera)
1008 {
1009 // FIXME: MHCollectionArea Trigger Area!
1010 if (header.IsDataRun())
1011 tasks.AddToList(&fillh3); // fill histogram task
1012 tasks.AddToList(&filltp); // fill histogram task
1013 }
1014 if (header.IsDataRun())
1015 {
1016 tasks.AddToList(&fillew); // fill histogram task
1017 tasks.AddToList(&filled); // fill histogram task
1018 }
1019 if (!header.IsPedestalRun())
1020 {
1021 tasks.AddToList(&fillx0a); // fill histogram task
1022 tasks.AddToList(&fillx0c); // fill histogram task
1023 }
1024 if (header.IsDataRun())
1025 {
1026 tasks.AddToList(&clean); // Cleaning for standard analysis
1027 tasks.AddToList(&hcalc); // Hillas parameter calculation
1028 tasks.AddToList(&cut);
1029 tasks.AddToList(&fillx0d); // fill histogram task
1030 tasks.AddToList(&fillx1); // fill histogram task
1031 //tasks.AddToList(&fillx2);
1032 tasks.AddToList(&fillx3); // fill histogram task
1033 tasks.AddToList(&fillx4); // fill histogram task
1034 if (!HasNullOut())
1035 tasks.AddToList(&write4a);
1036 //tasks.AddToList(&fillx5);
1037
1038 tasks.AddToList(&fillh4); // fill histogram task
1039 tasks.AddToList(&fillth); // fill histogram task
1040 tasks.AddToList(&fillca); // fill histogram task
1041 }
1042
1043
1044 // --------------------------------------------------------------------------------
1045 // --------------------------------------------------------------------------------
1046 // Event loop and processing display
1047 // --------------------------------------------------------------------------------
1048 // --------------------------------------------------------------------------------
1049
1050 // --------------------------------------------------------------------------------
1051 // Creation of event loop
1052 // --------------------------------------------------------------------------------
1053 tasks.SetAccelerator(MTask::kAccDontReset|MTask::kAccDontTime);
1054
1055 // Create and setup the eventloop
1056 MEvtLoop evtloop(fName);
1057 evtloop.SetParList(&plist);
1058 evtloop.SetDisplay(fDisplay);
1059 evtloop.SetLogStream(&gLog);
1060 if (!SetupEnv(evtloop))
1061 return kFALSE;
1062
1063 // FIXME: If pedestal file given use the values from this file
1064
1065
1066 // --------------------------------------------------------------------------------
1067 // Preparation of the graphicalk display which is shown while processing
1068 // --------------------------------------------------------------------------------
1069
1070 MGeomCam *cam = static_cast<MGeomCam*>(env2.GetCont());
1071
1072 write3af.SetGeometry(cam->ClassName());
1073
1074 MBinning *binsd = (MBinning*) plist.FindObject("BinningDistC");
1075 MBinning *binsd0 = (MBinning*) plist.FindObject("BinningDist");
1076 MBinning *binstr = (MBinning*) plist.FindObject("BinningTrigPos");
1077 if (binstr->IsDefault())
1078 binstr->SetEdgesLin(150, -header.GetFreqSampling()/1000.*shape.GetWidth(),
1079 header.GetFreqSampling()/1000.*(header.GetNumSamples()+shape.GetWidth()));
1080
1081 if (binsd->IsDefault() && cam)
1082 binsd->SetEdgesLin(100, 0, cam->GetMaxRadius()*cam->GetConvMm2Deg());
1083
1084 if (binsd0->IsDefault() && cam)
1085 binsd0->SetEdgesLin(100, 0, cam->GetMaxRadius()*cam->GetConvMm2Deg());
1086
1087
1088 header.Print();
1089
1090 // FIXME: Display acceptance curves!
1091
1092 if (fDisplay)
1093 {
1094 TCanvas &c = fDisplay->AddTab("Info");
1095 c.Divide(2,2);
1096
1097 c.cd(1);
1098 gPad->SetBorderMode(0);
1099 gPad->SetFrameBorderMode(0);
1100 gPad->SetGridx();
1101 gPad->SetGridy();
1102 gROOT->SetSelectedPad(0);
1103 shape.DrawClone()->SetBit(kCanDelete);
1104
1105 if (env0.GetCont() && (header.IsDataRun() || header.IsPointRun()))
1106 {
1107 c.cd(3);
1108 gPad->SetBorderMode(0);
1109 gPad->SetFrameBorderMode(0);
1110 gPad->SetGridx();
1111 gPad->SetGridy();
1112 gROOT->SetSelectedPad(0);
1113 static_cast<MOptics*>(env0.GetCont())->DrawClone("line")->SetBit(kCanDelete);
1114 }
1115
1116 if (fCamera)
1117 {
1118 if (env1.GetCont())
1119 {
1120 c.cd(2);
1121 gPad->SetBorderMode(0);
1122 gPad->SetFrameBorderMode(0);
1123 gPad->SetGridx();
1124 gPad->SetGridy();
1125 gROOT->SetSelectedPad(0);
1126 MHCamera *hcam = new MHCamera(static_cast<MGeomCam&>(*env1.GetCont()));
1127 hcam->SetStats(kFALSE);
1128 hcam->SetBit(MHCamera::kNoLegend);
1129 hcam->SetBit(kCanDelete);
1130 hcam->Draw();
1131 }
1132
1133 if (cam)
1134 {
1135 c.cd(4);
1136 gPad->SetBorderMode(0);
1137 gPad->SetFrameBorderMode(0);
1138 gPad->SetGridx();
1139 gPad->SetGridy();
1140 gROOT->SetSelectedPad(0);
1141 MHCamera *hcam = new MHCamera(*cam);
1142 hcam->SetStats(kFALSE);
1143 hcam->SetBit(MHCamera::kNoLegend);
1144 hcam->SetBit(kCanDelete);
1145 hcam->Draw();
1146 }
1147 }
1148
1149 TCanvas &d = fDisplay->AddTab("Info2");
1150 d.Divide(2,3);
1151
1152 d.cd(1);
1153 gPad->SetBorderMode(0);
1154 gPad->SetFrameBorderMode(0);
1155 gPad->SetGridx();
1156 gPad->SetGridy();
1157 gROOT->SetSelectedPad(0);
1158 splinepde.DrawClone()->SetBit(kCanDelete);
1159
1160 d.cd(3);
1161 gPad->SetBorderMode(0);
1162 gPad->SetFrameBorderMode(0);
1163 gPad->SetGridx();
1164 gPad->SetGridy();
1165 gROOT->SetSelectedPad(0);
1166 splinecones2.DrawClone()->SetBit(kCanDelete);
1167
1168 d.cd(5);
1169 gPad->SetBorderMode(0);
1170 gPad->SetFrameBorderMode(0);
1171 gPad->SetGridx();
1172 gPad->SetGridy();
1173 gROOT->SetSelectedPad(0);
1174 splinemirror.DrawClone()->SetBit(kCanDelete);
1175
1176 d.cd(2);
1177 gPad->SetBorderMode(0);
1178 gPad->SetFrameBorderMode(0);
1179 gPad->SetGridx();
1180 gPad->SetGridy();
1181 gROOT->SetSelectedPad(0);
1182 splinecones.DrawClone()->SetBit(kCanDelete);
1183 //splinecones2.DrawClone("same")->SetBit(kCanDelete);
1184
1185 d.cd(6);
1186 gPad->SetBorderMode(0);
1187 gPad->SetFrameBorderMode(0);
1188 gPad->SetGridx();
1189 gPad->SetGridy();
1190 gROOT->SetSelectedPad(0);
1191 MParSpline *sall = (MParSpline*)splinepde.DrawClone();
1192 //all->SetTitle("Combined acceptance");
1193 sall->SetBit(kCanDelete);
1194 if (splinemirror.GetSpline())
1195 sall->Multiply(*splinemirror.GetSpline());
1196 if (splinecones2.GetSpline())
1197 sall->Multiply(*splinecones2.GetSpline());
1198 }
1199
1200 // --------------------------------------------------------------------------------
1201 // Perform event loop
1202 // --------------------------------------------------------------------------------
1203
1204 // Execute first analysis
1205 if (!evtloop.Eventloop(fMaxEvents))
1206 {
1207 gLog << err << GetDescriptor() << ": Failed." << endl;
1208 return kFALSE;
1209 }
1210
1211 //-------------------------------------------
1212 // FIXME: Display Spline in tab
1213 // FIXME: Display Reflector in tab
1214 // FIXME: Display Routing(?) in tab
1215 // FIXME: Display Camera(s) in tab
1216 //-------------------------------------------
1217
1218 if (!WriteResult(plist, seq, header.GetRunNumber()))
1219 return kFALSE;
1220
1221 *fLog << all << GetDescriptor() << ": Done." << endl << endl << endl;;
1222
1223 return kTRUE;
1224}
Note: See TracBrowser for help on using the repository browser.