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

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