source: branches/Corsika7405Compatibility/mjobs/MJSimulation.cc

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