source: branches/Mars_MC/mjobs/MJSimulation.cc@ 17011

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