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

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