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