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

Last change on this file since 17559 was 17559, checked in by aparavac, 11 years ago
Added compatibilty for gcc 4.7 and c11 std
File size: 35.9 KB
Line 
1/* ======================================================================== *\
2!
3! *
4! * This file is part of MARS, the MAGIC Analysis and Reconstruction
5! * Software. It is distributed to you in the hope that it can be a useful
6! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
7! * It is distributed WITHOUT ANY WARRANTY.
8! *
9! * Permission to use, copy, modify and distribute this software and its
10! * documentation for any purpose is hereby granted without fee,
11! * provided that the above copyright notice appear in all copies and
12! * that both that copyright notice and this permission notice appear
13! * in supporting documentation. It is provided "as is" without express
14! * or implied warranty.
15! *
16!
17!
18! Author(s): Thomas Bretz, 1/2009 <mailto:tbretz@astro.uni-wuerzburg.de>
19!
20! Copyright: MAGIC Software Development, 2000-2009
21!
22!
23\* ======================================================================== */
24
25/////////////////////////////////////////////////////////////////////////////
26//
27// MJSimulation
28//
29//
30// Force reading a corsika file even if the footer (RUNE-section) is missing
31// by setting fForceMode to kTRUE or from the resource file by
32//
33// ForceMode: Yes
34//
35//
36// To switch off the simulation of the camera electronics, use:
37//
38// Camera: Off
39//
40// Note, that the border of the camera and the propertied of the cones
41// are still simulated (simply because it is fast). Furthermore, this
42// switches off the trigger for the output, i.e. all data which deposits
43// at least one photon in at least one pixel is written to the output file.
44//
45//
46// In case of a pedestal or calibration run the artificial trigger can
47// be "switched off" and the cosmics trrigger "switched on" by setting
48// fForceTrigger to kTRUE or from the resource file by
49//
50// ForceTrigger: Yes
51//
52//
53/////////////////////////////////////////////////////////////////////////////
54#include "MJSimulation.h"
55
56#include <TEnv.h>
57#include <TCanvas.h>
58
59// Core
60#include "MLog.h"
61#include "MLogManip.h"
62
63#include "MArgs.h"
64#include "MDirIter.h"
65#include "MParList.h"
66#include "MTaskList.h"
67#include "MEvtLoop.h"
68
69#include "MStatusDisplay.h"
70
71// Tasks
72#include "MCorsikaRead.h"
73#include "MContinue.h"
74#include "MFillH.h"
75#include "MGeomApply.h"
76#include "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
231void MJSimulation::SetupCommonFileStructure(MWriteRootFile &write) const
232{
233 // Common run headers
234 write.AddContainer("MMcCorsikaRunHeader", "RunHeaders", kFALSE);
235 write.AddContainer("MCorsikaRunHeader", "RunHeaders", kFALSE);
236 write.AddContainer("MRawRunHeader", "RunHeaders");
237 write.AddContainer("MGeomCam", "RunHeaders");
238 write.AddContainer("MMcRunHeader", "RunHeaders");
239
240 // Common events
241 write.AddContainer("MCorsikaEvtHeader", "Events", kFALSE);
242 write.AddContainer("MRawEvtHeader", "Events");
243 write.AddContainer("MMcEvt", "Events");
244 write.AddContainer("IncidentAngle", "Events", kFALSE);
245}
246//FIXME Etienne. I'm doing the same for fits and root. I probably should adapt somehow
247void MJSimulation::SetupCommonFileStructure(MWriteFitsFile& write) const
248{
249 // Common run headers
250 write.AddContainer("MMcCorsikaRunHeader", "RunHeaders", kFALSE);
251 write.AddContainer("MCorsikaRunHeader", "RunHeaders", kFALSE);
252 write.AddContainer("MRawRunHeader", "RunHeaders");
253 write.AddContainer("MGeomCam", "RunHeaders");
254 write.AddContainer("MMcRunHeader", "RunHeaders");
255
256 // Common events
257 write.AddContainer("MCorsikaEvtHeader", "Events", kFALSE);
258 write.AddContainer("MRawEvtHeader", "Events");
259 write.AddContainer("MMcEvt", "Events");
260 write.AddContainer("IncidentAngle", "Events", kFALSE);
261}
262
263void MJSimulation::SetupHeaderKeys(MWriteFitsFile &write,MRawRunHeader &header) const
264{
265 const MTime now(-1);
266 write.SetHeaderKey("ISMC",true,"Bool if File is Montecarlo File");
267 write.SetHeaderKey("TELESCOP", "FACT", "");
268 write.SetHeaderKey("PACKAGE", "MARS Cheobs", "");
269 write.SetHeaderKey("VERSION", "1.0", "");
270 write.SetHeaderKey("CREATOR", "Ceres", "");
271 write.SetHeaderKey("EXTREL", 1., "");
272 write.SetHeaderKey("COMPILED", __DATE__" " __TIME__, "");
273 write.SetHeaderKey("REVISION", "0", "");
274 write.SetHeaderKey("ORIGIN", "FACT", "");
275 write.SetHeaderKey("DATE", now.GetStringFmt("%Y-%m-%dT%H:%M:%S").Data(), "");
276 write.SetHeaderKey("NIGHT", now.GetNightAsInt(), "");
277 write.SetHeaderKey("TIMESYS", "UTC", "");
278 write.SetHeaderKey("TIMEUNIT", "d", "");
279 write.SetHeaderKey("MJDREF", 40587, "");
280 //write.SetHeaderKey("BLDVER", 1, "");
281 write.SetHeaderKey("RUNID", header.GetRunNumber(), "");
282 write.SetHeaderKey("NBOARD", 40, "");
283 write.SetHeaderKey("NPIX", header.GetNumPixel(), "");
284 write.SetHeaderKey("NROI", header.GetNumSamplesHiGain(), "");
285 write.SetHeaderKey("NROITM", 0, "");
286 write.SetHeaderKey("TMSHIFT", 0, "");
287 write.SetHeaderKey("CAMERA", "MGeomCamFACT", "");
288 write.SetHeaderKey("DAQ", "DRS4", "");
289
290 // FTemme: ADCRANGE and ADC have to be calculated, using the values for
291 // the fadctype.
292// write.SetHeaderKey("ADCRANGE", 2000, "Dynamic range in mV");
293// write.SetHeaderKey("ADC", 12, "Resolution in bits");
294
295 switch(header.GetRunType())
296 {
297 case MRawRunHeader::kRTData|MRawRunHeader::kRTMonteCarlo:
298 write.SetHeaderKey("RUNTYPE", "data", "");
299 break;
300 case MRawRunHeader::kRTPedestal|MRawRunHeader::kRTMonteCarlo:
301 write.SetHeaderKey("RUNTYPE", "pedestal", "");
302 break;
303 case MRawRunHeader::kRTCalibration|MRawRunHeader::kRTMonteCarlo:
304 write.SetHeaderKey("RUNTYPE", "calibration", "");
305 break;
306 }
307// write.SetHeaderKey("ID", 777, "Board 0: Board ID");
308// write.SetHeaderKey("FMVER", 532, "Board 0: Firmware Version");
309// write.SetHeaderKey("DNA", "0", "");
310// write.SetHeaderKey("BOARD", 0, "");
311// write.SetHeaderKey("PRESC", 40, "");
312// write.SetHeaderKey("PHASE", 0, "");
313// write.SetHeaderKey("DAC0", 26500, "");
314// write.SetHeaderKey("DAC1", 0, "");
315// write.SetHeaderKey("DAC2", 0, "");
316// write.SetHeaderKey("DAC3", 0, "");
317// write.SetHeaderKey("DAC4", 28800, "");
318// write.SetHeaderKey("DAC5", 28800, "");
319// write.SetHeaderKey("DAC6", 28800, "");
320// write.SetHeaderKey("DAC7", 28800, "");
321 write.SetHeaderKey("REFCLK", header.GetFreqSampling(), "");
322 write.SetHeaderKey("DRSCALIB", false, "");
323// write.SetHeaderKey("TSTARTI", 0, "");
324// write.SetHeaderKey("TSTARTF", 0., "");
325// write.SetHeaderKey("TSTOPI", 0, "");
326// write.SetHeaderKey("TSTOPF", 0., "");
327// write.SetHeaderKey("DATE-OBS", "1970-01-01T00:00:00", "");
328// write.SetHeaderKey("DATE-END", "1970-01-01T00:00:00", "");
329// write.SetHeaderKey("NTRG", 0, "");
330// write.SetHeaderKey("NTRGPED", 0, "");
331// write.SetHeaderKey("NTRGLPE", 0, "");
332// write.SetHeaderKey("NTRGTIM", 0, "");
333// write.SetHeaderKey("NTRGLPI", 0, "");
334// write.SetHeaderKey("NTRGEXT1", 0, "");
335// write.SetHeaderKey("NTRGEXT2", 0, "");
336// write.SetHeaderKey("NTRGMISC", 0, "");
337}
338
339void MJSimulation::SetupVetoColumns(MWriteFitsFile &write) const
340{
341 write.VetoColumn("MParameterD.fVal");
342 write.VetoColumn("MRawEvtData.fLoGainPixId");
343 write.VetoColumn("MRawEvtData.fLoGainFadcSamples");
344 write.VetoColumn("MRawEvtData.fABFlags");
345 write.VetoColumn("MRawEvtHeader.fNumTrigLvl2");
346 //write.VetoColumn("MRawEvtHeader.fTrigPattern");
347 write.VetoColumn("MRawEvtHeader.fNumLoGainOn");
348}
349
350Bool_t MJSimulation::Process(const MArgs &args, const MSequence &seq)
351{
352 /*
353 if (!fSequence.IsValid())
354 {
355 *fLog << err << "ERROR - Sequence invalid!" << endl;
356 return kFALSE;
357 }
358 */
359
360// if (!HasWritePermission(CombinePath(fPathOut, GetOutFile(seq))))
361// return kFALSE;
362
363 *fLog << inf;
364 fLog->Separator(GetDescriptor());
365
366 if (!CheckEnv())
367 return kFALSE;
368
369 if (seq.IsValid())
370 *fLog << fSequence.GetFileName() << endl;
371 else
372 *fLog << "Input: " << args.GetNumArguments() << "-files" << endl;
373 *fLog << endl;
374
375 MDirIter iter;
376 if (seq.IsValid() && seq.GetRuns(iter, MSequence::kCorsika)<=0)
377 {
378 *fLog << err << "ERROR - Sequence valid but without files." << endl;
379 return kFALSE;
380 }
381
382 // --------------------------------------------------------------------------------
383
384 // Setup Parlist
385 MParList plist;
386 plist.AddToList(this); // take care of fDisplay!
387
388 // setup TaskList
389 MTaskList tasks;
390 plist.AddToList(&tasks);
391
392 // --------------------------------------------------------------------------------
393
394 // FIXME: Allow empty GeomCones!
395 MParEnv env0("Reflector");
396 MParEnv env1("GeomCones");
397 MParEnv env2("MGeomCam");
398 env0.SetClassName("MReflector");
399 env1.SetClassName("MGeomCam");
400 env2.SetClassName("MGeomCam");
401 plist.AddToList(&env0);
402 plist.AddToList(&env1);
403 plist.AddToList(&env2);
404
405 plist.FindCreateObj("MPedPhotCam", "MPedPhotFromExtractorRndm");
406
407 MParSpline shape("PulseShape");
408 plist.AddToList(&shape);
409
410 // *** FIXME *** FIXME *** FIXME ***
411 plist.FindCreateObj("MMcRunHeader");
412
413 MRawRunHeader header;
414 header.SetValidMagicNumber();
415 //header.InitFadcType(3);
416
417 switch (fOperationMode)
418 {
419 case kModeData:
420 header.SetRunType(MRawRunHeader::kRTMonteCarlo|MRawRunHeader::kRTData);
421 header.SetRunInfo(0, fRunNumber<0 ? 3 : fRunNumber);
422 break;
423
424 case kModePed:
425 header.SetRunType(MRawRunHeader::kRTMonteCarlo|MRawRunHeader::kRTPedestal);
426 header.SetSourceInfo("Pedestal");
427 header.SetRunInfo(0, fRunNumber<0 ? 1 : fRunNumber);
428 break;
429
430 case kModeCal:
431 header.SetRunType(MRawRunHeader::kRTMonteCarlo|MRawRunHeader::kRTCalibration);
432 header.SetSourceInfo("Calibration");
433 header.SetRunInfo(0, fRunNumber<0 ? 2 : fRunNumber);
434 break;
435
436 case kModePointRun:
437 header.SetRunType(MRawRunHeader::kRTMonteCarlo|MRawRunHeader::kRTPointRun);
438 header.SetRunInfo(0, fRunNumber<0 ? 0 : fRunNumber);
439 break;
440 }
441
442 // FIXME: Move to MSimPointingPos, MSimCalibrationSignal
443 // Can we use this as input for MSimPointingPos?
444 header.SetObservation("On", "MonteCarlo");
445 plist.AddToList(&header);
446 // ++++++++ FIXME FIXME FIXME +++++++++++++
447
448 /*
449 MPedestalCam pedcam;
450 pedcam.Init(geomcam.GetNumPixels());
451 for (UInt_t i=0; i<geomcam.GetNumPixels(); i++)
452 pedcam[i].Set(128./header.GetScale(), 22.5/header.GetScale());
453 plist.AddToList(&pedcam);
454 */
455
456 // -------------------------------------------------------------------
457
458 MCorsikaRead read;
459 read.SetForceMode(fForceMode);
460
461 if (!seq.IsValid())
462 {
463 for (int i=0; i<args.GetNumArguments(); i++)
464 read.AddFile(args.GetArgumentStr(i));
465 }
466 else
467 read.AddFiles(iter);
468
469 MContinue precut("", "PreCut");
470 precut.IsInverted();
471 precut.SetAllowEmpty();
472
473 MSimMMCS simmmcs;
474
475 MParSpline splinepde("PhotonDetectionEfficiency");
476 MParSpline splinemirror("MirrorReflectivity");
477 MParSpline splinecones("ConesAngularAcceptance");
478 MParSpline splinecones2("ConesTransmission");
479 plist.AddToList(&splinepde);
480 plist.AddToList(&splinemirror);
481 plist.AddToList(&splinecones);
482 plist.AddToList(&splinecones2);
483
484 const TString sin2 = "sin(MCorsikaEvtHeader.fZd)*sin(MCorsikaRunHeader.fZdMin*TMath::DegToRad())";
485 const TString cos2 = "cos(MCorsikaEvtHeader.fZd)*cos(MCorsikaRunHeader.fZdMin*TMath::DegToRad())";
486 const TString cos = "cos(MCorsikaEvtHeader.fAz-MCorsikaRunHeader.fAzMin*TMath::DegToRad())";
487
488 const TString form = "acos("+sin2+"*"+cos+"+"+cos2+")*TMath::RadToDeg()";
489
490 MParameterCalc calcangle(form, "CalcIncidentAngle");
491 calcangle.SetNameParameter("IncidentAngle");
492
493 MSimAtmosphere simatm;
494 MSimAbsorption absapd("SimPhotonDetectionEfficiency");
495 MSimAbsorption absmir("SimMirrorReflectivity");
496 MSimAbsorption cones("SimConesAngularAcceptance");
497 MSimAbsorption cones2("SimConesTransmission");
498 absapd.SetParName("PhotonDetectionEfficiency");
499 absmir.SetParName("MirrorReflectivity");
500 cones.SetParName("ConesAngularAcceptance");
501 cones.SetUseTheta();
502 cones2.SetParName("ConesTransmission");
503
504 MSimPointingPos pointing;
505
506 MSimReflector reflect;
507 reflect.SetNameGeomCam("GeomCones");
508 reflect.SetNameReflector("Reflector");
509// MSimStarField stars;
510
511 MContinue cont1("MPhotonEvent.GetNumPhotons<1", "ContEmpty1");
512 MContinue cont2("MPhotonEvent.GetNumPhotons<1", "ContEmpty2");
513 MContinue cont3("MPhotonEvent.GetNumPhotons<2", "ContEmpty3");
514 cont1.SetAllowEmpty();
515 cont2.SetAllowEmpty();
516 cont3.SetAllowEmpty();
517
518 // -------------------------------------------------------------------
519
520 MBinning binse( 120, 1, 1000000, "BinningEnergy", "log");
521 MBinning binsth( 60, 0.9, 900000, "BinningThreshold", "log");
522 MBinning binsee( 36, 0.9, 900000, "BinningEnergyEst", "log");
523 MBinning binss( 100, 1, 10000000, "BinningSize", "log");
524// MBinning binsi( 100, -500, 500, "BinningImpact");
525 MBinning binsi( 32, 0, 800, "BinningImpact");
526 MBinning binsh( 150, 0, 50, "BinningHeight");
527 MBinning binsaz(720, -360, 360, "BinningAz");
528 MBinning binszd( 70, 0, 70, "BinningZd");
529 MBinning binsvc( 35, 0, 7, "BinningViewCone");
530 MBinning binsel(150, 0, 50, "BinningTotLength");
531 MBinning binsew(150, 0, 15, "BinningMedLength");
532 MBinning binsd("BinningDistC");
533 MBinning binsd0("BinningDist");
534 MBinning binstr("BinningTrigPos");
535
536 plist.AddToList(&binsee);
537 plist.AddToList(&binse);
538 plist.AddToList(&binsth);
539 plist.AddToList(&binss);
540 plist.AddToList(&binsi);
541 plist.AddToList(&binsh);
542 plist.AddToList(&binszd);
543 plist.AddToList(&binsaz);
544 plist.AddToList(&binsvc);
545 plist.AddToList(&binsel);
546 plist.AddToList(&binsew);
547 plist.AddToList(&binstr);
548 plist.AddToList(&binsd);
549 plist.AddToList(&binsd0);
550
551 MHn mhn1, mhn2, mhn3, mhn4;
552 SetupHist(mhn1);
553 SetupHist(mhn2);
554 SetupHist(mhn3);
555 SetupHist(mhn4);
556
557 MH3 mhtp("TriggerPos.fVal-IntendedPulsePos.fVal-PulseShape.GetWidth");
558 mhtp.SetName("TrigPos");
559 mhtp.SetTitle("Trigger position w.r.t. the first photon hitting a detector");
560
561 MH3 mhew("MPhotonStatistics.fLength");
562 mhew.SetName("TotLength");
563 mhew.SetTitle("Time between first and last photon hitting a detector;L [ns]");
564
565 MH3 mhed("MPhotonStatistics.fTimeMedDev");
566 mhed.SetName("MedLength");
567 mhed.SetTitle("Median deviation (1\\sigma);L [ns]");
568
569 MFillH fillh1(&mhn1, "", "FillCorsika");
570 MFillH fillh2(&mhn2, "", "FillH2");
571 MFillH fillh3(&mhn3, "", "FillH3");
572 MFillH fillh4(&mhn4, "", "FillH4");
573 MFillH filltp(&mhtp, "", "FillTriggerPos");
574 MFillH fillew(&mhew, "", "FillEvtWidth");
575 MFillH filled(&mhed, "", "FillMedDev");
576 fillh1.SetNameTab("Corsika", "Distribution as simulated by Corsika");
577 fillh2.SetNameTab("Detectable", "Distribution of events hit the detector");
578 fillh3.SetNameTab("Triggered", "Distribution of triggered events");
579 fillh4.SetNameTab("Cleaned", "Distribution after cleaning and cuts");
580 filltp.SetNameTab("TrigPos", "TriggerPosition w.r.t the first photon");
581 fillew.SetNameTab("EvtWidth", "Time between first and last photon hitting a detector");
582 filled.SetNameTab("MedDev", "Median deviation of arrival time of photons hitting a detector");
583
584 MHPhotonEvent planeG(1, "HPhotonEventGround"); // Get from MaxImpact
585 MHPhotonEvent plane0(2, "HMirrorPlane0"); // Get from MReflector
586 //MHPhotonEvent plane1(2, "HMirrorPlane1");
587 MHPhotonEvent plane2(2, "HMirrorPlane2");
588 MHPhotonEvent plane3(2, "HMirrorPlane3");
589 MHPhotonEvent plane4(2, "HMirrorPlane4");
590 MHPhotonEvent planeF1(6, "HPhotonEventCamera"); // Get from MGeomCam
591 MHPhotonEvent planeF2(header.IsPointRun()?4:6, "HPhotonEventCones"); // Get from MGeomCam
592
593 plist.AddToList(&planeG);
594 plist.AddToList(&plane0);
595 plist.AddToList(&plane2);
596 plist.AddToList(&plane3);
597 plist.AddToList(&plane4);
598 plist.AddToList(&planeF1);
599 plist.AddToList(&planeF2);;
600
601 //MHPSF psf;
602
603 MFillH fillG(&planeG, "MPhotonEvent", "FillGround");
604 MFillH fill0(&plane0, "MirrorPlane0", "FillReflector");
605 //MFillH fill1(&plane1, "MirrorPlane1", "FillCamShadow");
606 MFillH fill2(&plane2, "MirrorPlane2", "FillCandidates");
607 MFillH fill3(&plane3, "MirrorPlane3", "FillReflected");
608 MFillH fill4(&plane4, "MirrorPlane4", "FillFocal");
609 MFillH fillF1(&planeF1, "MPhotonEvent", "FillCamera");
610 MFillH fillF2(&planeF2, "MPhotonEvent", "FillCones");
611 //MFillH fillP(&psf, "MPhotonEvent", "FillPSF");
612
613 fillG.SetNameTab("Ground", "Photon distribution at ground");
614 fill0.SetNameTab("Reflector", "Photon distribution at reflector plane w/o camera shadow");
615 //fill1.SetNameTab("CamShadow", "Photon distribution at reflector plane w/ camera shadow");
616 fill2.SetNameTab("Candidates", "*Can hit* photon distribution at reflector plane w/ camera shadow");
617 fill3.SetNameTab("Reflected", "Photon distribution after reflector projected to reflector plane");
618 fill4.SetNameTab("Focal", "Photon distribution at focal plane");
619 fillF1.SetNameTab("Camera", "Photon distribution which hit the detector");
620 fillF2.SetNameTab("Cones", "Photon distribution after cones");
621 //fillP.SetNameTab("PSF", "Photon distribution after cones for all mirrors");
622
623 // -------------------------------------------------------------------
624
625 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());
626
627 const TString rule1(Form(fmt, 'R'));
628 const TString rule2(Form(fmt, 'Y'));
629 const TString rule4(Form(fmt, 'I'));
630 TString rule3(Form(fmt, header.GetRunTypeChar()));
631
632 MWriteRootFile write4a( 2, rule4, fOverwrite?"RECREATE":"NEW", "Star file");
633 MWriteRootFile write4b( 2, rule4, fOverwrite?"RECREATE":"NEW", "Star file");
634 MWriteRootFile write3b( 2, rule3, fOverwrite?"RECREATE":"NEW", "Camera file");
635 MWriteRootFile write2a( 2, rule2, fOverwrite?"RECREATE":"NEW", "Signal file");
636 MWriteRootFile write2b( 2, rule2, fOverwrite?"RECREATE":"NEW", "Signal file");
637 MWriteRootFile write1a( 2, rule1, fOverwrite?"RECREATE":"NEW", "Reflector file");
638 MWriteRootFile write1b( 2, rule1, fOverwrite?"RECREATE":"NEW", "Reflector file");
639
640 if (fWriteFitsFile)
641 rule3.ReplaceAll("$1.root/", "$1.fits/");
642
643 MWriteFitsFile write3af( 2, rule3, fOverwrite?"RECREATE":"NEW", "Camera file");
644 MWriteRootFile write3ar( 2, rule3, fOverwrite?"RECREATE":"NEW", "Camera file");
645
646 MTask &write3a = fWriteFitsFile ? static_cast<MTask&>(write3af) : static_cast<MTask&>(write3ar);
647
648 //SetupHeaderKeys(write3af,header);
649 SetupVetoColumns(write3af);
650
651 write3af.SetBytesPerSample("Data", 2);
652
653 write1a.SetName("WriteRefData");
654 write1b.SetName("WriteRefMC");
655 write2a.SetName("WriteSigData");
656 write2b.SetName("WriteSigMC");
657 write3a.SetName("WriteCamData");
658 write3b.SetName("WriteCamMC");
659 write4a.SetName("WriteImgData");
660 write4b.SetName("WriteImgMC");
661
662 SetupCommonFileStructure(write1a);
663 SetupCommonFileStructure(write2a);
664 SetupCommonFileStructure(write3ar);
665 SetupCommonFileStructure(write3af);
666 SetupCommonFileStructure(write4a);
667
668 // R: Dedicated file structure
669 write1a.AddContainer("MPhotonEvent", "Events");
670
671 // Y: Dedicated file structure
672 write2a.AddContainer("MPedPhotFromExtractorRndm", "RunHeaders"); // FIXME: Needed for the signal files to be display in MARS
673 write2a.AddContainer("MSignalCam", "Events");
674
675 // D: Dedicated file structure
676 write3af.AddContainer("ElectronicNoise", "RunHeaders");
677 write3af.AddContainer("IntendedPulsePos", "RunHeaders");
678 write3af.AddContainer("MRawEvtData", "Events");
679
680 if (fOperationMode == kModeData){
681 write3af.AddContainer("MMcEvtBasic", "Events");
682 }
683
684 write3af.AddContainer("MPhotonEvent","Events");
685
686 write3ar.AddContainer("ElectronicNoise", "RunHeaders");
687 write3ar.AddContainer("IntendedPulsePos", "RunHeaders");
688 write3ar.AddContainer("MRawEvtData", "Events");
689 // It doesn't make much sene to write this information
690 // to the file because the units are internal units not
691 // related to the output samples
692 // if (header.IsDataRun() || fForceTrigger)
693 // write3a.AddContainer("TriggerPos", "Events");
694
695 // I: Dedicated file structure
696 write4a.AddContainer("MHillas", "Events");
697 write4a.AddContainer("MHillasSrc", "Events");
698 write4a.AddContainer("MImagePar", "Events");
699 write4a.AddContainer("MNewImagePar", "Events");
700
701 // Basic MC data
702 write1b.AddContainer("MMcEvtBasic", "OriginalMC");
703 write2b.AddContainer("MMcEvtBasic", "OriginalMC");
704 write3b.AddContainer("MMcEvtBasic", "OriginalMC");
705 write4b.AddContainer("MMcEvtBasic", "OriginalMC");
706
707 // -------------------------------------------------------------------
708
709 MGeomApply apply;
710
711 MSimPSF simpsf;
712
713 MSimGeomCam simgeom;
714 simgeom.SetNameGeomCam("GeomCones");
715
716 MSimCalibrationSignal simcal;
717 simcal.SetNameGeomCam("GeomCones");
718
719 // Sky Quality Meter: 21.82M = 2.02e-4cd/m^2
720 // 1cd = 12.566 lm / 4pi sr
721
722 // FIXME: Simulate photons before cones and QE!
723
724 MSimRandomPhotons simnsb;
725 simnsb.SetFreq(5.8, 0.004); // ph/m^2/nm/sr/ns NSB, 4MHz dark count rate
726 simnsb.SetNameGeomCam("GeomCones");
727
728 // FIXME: How do we handle star-light? We need the rate but we also
729 // need to process it through the mirror!
730
731 MSimAPD simapd;
732 simapd.SetNameGeomCam("GeomCones");
733
734 MSimExcessNoise simexcnoise;
735 MSimBundlePhotons simsum;
736 MSimSignalCam simsignal;
737 MSimCamera simcam;
738 MSimTrigger simtrig;
739 MSimReadout simdaq;
740
741 MContinue conttrig("TriggerPos.fVal<0", "ContTrigger");
742
743 MParameterD pulpos("IntendedPulsePos");
744 // FIXME: Set a default which could be 1/3 of the digitization window
745 // pulpos.SetVal(40); // [ns]
746 plist.AddToList(&pulpos);
747
748 MParameterD trigger("TriggerPos");
749 trigger.SetVal(0);
750 plist.AddToList(&trigger);
751
752 // -------------------------------------------------------------------
753
754 // Remove isolated pixels
755 MImgCleanStd clean(0, 0);
756 //clean.SetCleanLvl0(0); // The level above which isolated pixels are kept
757 clean.SetCleanRings(0);
758 clean.SetMethod(MImgCleanStd::kAbsolute);
759
760 //clean.SetNamePedPhotCam("MPedPhotFromExtractorRndm");
761
762 MHillasCalc hcalc;
763 hcalc.Disable(MHillasCalc::kCalcConc);
764
765 MHCamEvent evt0a(/*10*/3, "Signal", "Average signal (absolute);;S [ph]");
766 MHCamEvent evt0c(/*10*/3, "MaxSignal", "Maximum signal (absolute);;S [ph]");
767 MHCamEvent evt0d(/*11*/8, "ArrTm", "Time after first photon;;T [ns]");
768 evt0a.SetErrorSpread(kFALSE);
769 evt0c.SetCollectMax();
770
771 MContinue cut("", "Cut");
772
773 MFillH fillx0a(&evt0a, "MSignalCam", "FillAvgSignal");
774 MFillH fillx0c(&evt0c, "MSignalCam", "FillMaxSignal");
775 MFillH fillx0d(&evt0d, "MSignalCam", "FillArrTm");
776 MFillH fillx1("MHHillas", "MHillas", "FillHillas");
777 MFillH fillx3("MHHillasSrc", "MHillasSrc", "FillHillasSrc");
778 MFillH fillx4("MHNewImagePar", "MNewImagePar", "FillNewImagePar");
779 MFillH fillth("MHThreshold", "", "FillThreshold");
780 MFillH fillca("MHCollectionArea", "", "FillTrigArea");
781
782 fillth.SetNameTab("Threshold");
783 fillca.SetNameTab("TrigArea");
784
785 // -------------------------------------------------------------------
786
787 if (header.IsDataRun())
788 {
789 tasks.AddToList(&read);
790 tasks.AddToList(&precut);
791 tasks.AddToList(&pointing);
792 tasks.AddToList(&simmmcs);
793 if (!fPathOut.IsNull() && !HasNullOut())
794 {
795 //tasks.AddToList(&write1b);
796 tasks.AddToList(&write2b);
797 if (fCamera)
798 tasks.AddToList(&write3b);
799 if (header.IsDataRun())
800 tasks.AddToList(&write4b);
801 }
802 // if (header.IsPointRun())
803 // tasks.AddToList(&stars);
804 if (1)
805 tasks.AddToList(&simatm); // Here because before fillh1
806 tasks.AddToList(&calcangle);
807 tasks.AddToList(&fillh1);
808 tasks.AddToList(&fillG);
809 if (!header.IsPointRun())
810 {
811 tasks.AddToList(&absapd);
812 tasks.AddToList(&absmir);
813 if (0)
814 tasks.AddToList(&simatm); // FASTER?
815 }
816 tasks.AddToList(&reflect);
817 if (!header.IsPointRun())
818 {
819 tasks.AddToList(&fill0);
820 //tasks.AddToList(&fill1);
821 tasks.AddToList(&fill2);
822 tasks.AddToList(&fill3);
823 tasks.AddToList(&fill4);
824 tasks.AddToList(&fillF1);
825 }
826 tasks.AddToList(&cones);
827 tasks.AddToList(&cones2);
828 //if (header.IsPointRun())
829 // tasks.AddToList(&fillP);
830 tasks.AddToList(&cont1);
831 }
832 // -------------------------------
833 tasks.AddToList(&apply);
834 if (header.IsDataRun())
835 {
836 tasks.AddToList(&simpsf);
837 tasks.AddToList(&simgeom);
838 tasks.AddToList(&cont2);
839 if (!header.IsPointRun())
840 {
841 tasks.AddToList(&fillF2);
842 tasks.AddToList(&fillh2);
843 }
844 tasks.AddToList(&cont3);
845 }
846 if (fCamera)
847 {
848 if (header.IsPedestalRun() || header.IsCalibrationRun())
849 tasks.AddToList(&simcal);
850 tasks.AddToList(&simnsb);
851 tasks.AddToList(&simapd);
852 tasks.AddToList(&simexcnoise);
853 }
854 tasks.AddToList(&simsum);
855 if (fCamera)
856 {
857 tasks.AddToList(&simcam);
858 if (header.IsDataRun() || fForceTrigger)
859 tasks.AddToList(&simtrig);
860 tasks.AddToList(&conttrig);
861 tasks.AddToList(&simdaq);
862 }
863 tasks.AddToList(&simsignal); // What do we do if signal<0?
864 if (!fPathOut.IsNull() && !HasNullOut())
865 {
866 //tasks.AddToList(&write1a);
867 if (!header.IsPedestalRun())
868 tasks.AddToList(&write2a);
869 if (fCamera)
870 tasks.AddToList(&write3a);
871 }
872 // -------------------------------
873 if (fCamera)
874 {
875 // FIXME: MHCollectionArea Trigger Area!
876 if (header.IsDataRun())
877 tasks.AddToList(&fillh3);
878 tasks.AddToList(&filltp);
879 }
880 if (header.IsDataRun())
881 {
882 tasks.AddToList(&fillew);
883 tasks.AddToList(&filled);
884 }
885 if (!header.IsPedestalRun())
886 {
887 tasks.AddToList(&fillx0a);
888 tasks.AddToList(&fillx0c);
889 }
890 if (header.IsDataRun())
891 {
892 tasks.AddToList(&clean);
893 tasks.AddToList(&hcalc);
894 tasks.AddToList(&cut);
895 tasks.AddToList(&fillx0d);
896 tasks.AddToList(&fillx1);
897 //tasks.AddToList(&fillx2);
898 tasks.AddToList(&fillx3);
899 tasks.AddToList(&fillx4);
900 if (!HasNullOut())
901 tasks.AddToList(&write4a);
902 //tasks.AddToList(&fillx5);
903
904 tasks.AddToList(&fillh4);
905 tasks.AddToList(&fillth);
906 tasks.AddToList(&fillca);
907 }
908 //-------------------------------------------
909
910 tasks.SetAccelerator(MTask::kAccDontReset|MTask::kAccDontTime);
911
912 // Create and setup the eventloop
913 MEvtLoop evtloop(fName);
914 evtloop.SetParList(&plist);
915 evtloop.SetDisplay(fDisplay);
916 evtloop.SetLogStream(&gLog);
917 if (!SetupEnv(evtloop))
918 return kFALSE;
919
920 // FIXME: If pedestal file given use the values from this file
921 //-------------------------------------------
922
923 MGeomCam *cam = static_cast<MGeomCam*>(env2.GetCont());
924
925 if (binstr.IsDefault())
926 binstr.SetEdgesLin(150, -shape.GetWidth(),
927 header.GetFreqSampling()/1000.*header.GetNumSamples()+shape.GetWidth());
928
929 if (binsd.IsDefault() && cam)
930 binsd.SetEdgesLin(100, 0, cam->GetMaxRadius()*cam->GetConvMm2Deg());
931
932 if (binsd0.IsDefault() && cam)
933 binsd0.SetEdgesLin(100, 0, cam->GetMaxRadius()*cam->GetConvMm2Deg());
934
935
936 header.Print();
937
938 // FIXME: Display acceptance curves!
939
940 if (fDisplay)
941 {
942 TCanvas &c = fDisplay->AddTab("Info");
943 c.Divide(2,2);
944
945 c.cd(1);
946 gPad->SetBorderMode(0);
947 gPad->SetFrameBorderMode(0);
948 gPad->SetGridx();
949 gPad->SetGridy();
950 gROOT->SetSelectedPad(0);
951 shape.DrawClone()->SetBit(kCanDelete);
952
953 if (env0.GetCont() && (header.IsDataRun() || header.IsPointRun()))
954 {
955 c.cd(3);
956 gPad->SetBorderMode(0);
957 gPad->SetFrameBorderMode(0);
958 gPad->SetGridx();
959 gPad->SetGridy();
960 gROOT->SetSelectedPad(0);
961 static_cast<MReflector*>(env0.GetCont())->DrawClone("line")->SetBit(kCanDelete);
962 }
963
964 if (fCamera)
965 {
966 if (env1.GetCont())
967 {
968 c.cd(2);
969 gPad->SetBorderMode(0);
970 gPad->SetFrameBorderMode(0);
971 gPad->SetGridx();
972 gPad->SetGridy();
973 gROOT->SetSelectedPad(0);
974 MHCamera *c = new MHCamera(static_cast<MGeomCam&>(*env1.GetCont()));
975 c->SetStats(kFALSE);
976 c->SetBit(MHCamera::kNoLegend);
977 c->SetBit(kCanDelete);
978 c->Draw();
979 }
980
981 if (cam)
982 {
983 c.cd(4);
984 gPad->SetBorderMode(0);
985 gPad->SetFrameBorderMode(0);
986 gPad->SetGridx();
987 gPad->SetGridy();
988 gROOT->SetSelectedPad(0);
989 MHCamera *c = new MHCamera(*cam);
990 c->SetStats(kFALSE);
991 c->SetBit(MHCamera::kNoLegend);
992 c->SetBit(kCanDelete);
993 c->Draw();
994 }
995 }
996
997 TCanvas &d = fDisplay->AddTab("Info2");
998 d.Divide(2,3);
999
1000 d.cd(1);
1001 gPad->SetBorderMode(0);
1002 gPad->SetFrameBorderMode(0);
1003 gPad->SetGridx();
1004 gPad->SetGridy();
1005 gROOT->SetSelectedPad(0);
1006 splinepde.DrawClone()->SetBit(kCanDelete);
1007
1008 d.cd(3);
1009 gPad->SetBorderMode(0);
1010 gPad->SetFrameBorderMode(0);
1011 gPad->SetGridx();
1012 gPad->SetGridy();
1013 gROOT->SetSelectedPad(0);
1014 splinecones2.DrawClone()->SetBit(kCanDelete);
1015
1016 d.cd(5);
1017 gPad->SetBorderMode(0);
1018 gPad->SetFrameBorderMode(0);
1019 gPad->SetGridx();
1020 gPad->SetGridy();
1021 gROOT->SetSelectedPad(0);
1022 splinemirror.DrawClone()->SetBit(kCanDelete);
1023
1024 d.cd(2);
1025 gPad->SetBorderMode(0);
1026 gPad->SetFrameBorderMode(0);
1027 gPad->SetGridx();
1028 gPad->SetGridy();
1029 gROOT->SetSelectedPad(0);
1030 splinecones.DrawClone()->SetBit(kCanDelete);
1031 //splinecones2.DrawClone("same")->SetBit(kCanDelete);
1032
1033 d.cd(6);
1034 gPad->SetBorderMode(0);
1035 gPad->SetFrameBorderMode(0);
1036 gPad->SetGridx();
1037 gPad->SetGridy();
1038 gROOT->SetSelectedPad(0);
1039 MParSpline *all = (MParSpline*)splinepde.DrawClone();
1040 //all->SetTitle("Combined acceptance");
1041 all->SetBit(kCanDelete);
1042 if (splinemirror.GetSpline())
1043 all->Multiply(*splinemirror.GetSpline());
1044 if (splinecones2.GetSpline())
1045 all->Multiply(*splinecones2.GetSpline());
1046 }
1047
1048 //-------------------------------------------
1049
1050 // Execute first analysis
1051 if (!evtloop.Eventloop(fMaxEvents))
1052 {
1053 gLog << err << GetDescriptor() << ": Failed." << endl;
1054 return kFALSE;
1055 }
1056
1057 //-------------------------------------------
1058 // FIXME: Display Spline in tab
1059 // FIXME: Display Reflector in tab
1060 // FIXME: Display Routing(?) in tab
1061 // FIXME: Display Camera(s) in tab
1062 //-------------------------------------------
1063
1064 if (!WriteResult(plist, seq, header.GetRunNumber()))
1065 return kFALSE;
1066
1067 *fLog << all << GetDescriptor() << ": Done." << endl << endl << endl;;
1068
1069 return kTRUE;
1070}
Note: See TracBrowser for help on using the repository browser.