source: trunk/MagicSoft/Mars/mjobs/MJSimulation.cc@ 9294

Last change on this file since 9294 was 9279, checked in by tbretz, 17 years ago
*** empty log message ***
File size: 19.4 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#include "MJSimulation.h"
31
32#include <TEnv.h>
33
34// Core
35#include "MLog.h"
36#include "MLogManip.h"
37
38#include "MArgs.h"
39//#include "MDirIter.h"
40#include "MParList.h"
41#include "MTaskList.h"
42#include "MEvtLoop.h"
43
44#include "MStatusDisplay.h"
45
46// Tasks
47#include "MCorsikaRead.h"
48#include "MContinue.h"
49#include "MFillH.h"
50#include "MGeomApply.h"
51#include "MHillasCalc.h"
52#include "MImgCleanStd.h"
53#include "MWriteRootFile.h"
54
55#include "MSimAbsorption.h"
56#include "MSimReflector.h"
57#include "MSimPointingPos.h"
58#include "MSimGeomCam.h"
59#include "MSimSignalCam.h"
60#include "MSimAPD.h"
61#include "MSimExcessNoise.h"
62#include "MSimCamera.h"
63#include "MSimTrigger.h"
64#include "MSimReadout.h"
65#include "MSimRandomPhotons.h"
66#include "MSimBundlePhotons.h"
67#include "MSimCalibrationSignal.h"
68
69// Histograms
70#include "MBinning.h"
71
72#include "MHn.h"
73#include "MHCamEvent.h"
74#include "MHPhotonEvent.h"
75
76// Container
77#include "MRawRunHeader.h"
78#include "MParameters.h"
79#include "MReflector.h"
80#include "MParEnv.h"
81#include "MPulseShape.h"
82
83#include "MPedestalCam.h"
84#include "MPedestalPix.h"
85
86ClassImp(MJSimulation);
87
88using namespace std;
89
90// --------------------------------------------------------------------------
91//
92// Default constructor.
93//
94// Sets fRuns to 0, fExtractor to NULL, fDataCheck to kFALSE
95//
96MJSimulation::MJSimulation(const char *name, const char *title) : fForceMode(kFALSE)
97{
98 fName = name ? name : "MJSimulation";
99 fTitle = title ? title : "Standard analysis and reconstruction";
100}
101
102Bool_t MJSimulation::CheckEnvLocal()
103{
104 fForceMode = GetEnv("ForceMode", fForceMode);
105
106 return kTRUE;
107}
108
109Bool_t MJSimulation::WriteResult()
110{
111 if (fPathOut.IsNull())
112 {
113 *fLog << inf << "No output path specified via SetPathOut - no output written." << endl;
114 return kTRUE;
115 }
116
117 TObjArray cont;
118 cont.Add(const_cast<TEnv*>(GetEnv()));
119 //cont.Add(const_cast<MSequence*>(&fSequence));
120
121 if (fDisplay)
122 {
123// TString title = "-- Reflector: ";
124// title += fSequence.GetSequence();
125// title += " --";
126// fDisplay->SetTitle(title, kFALSE);
127 fDisplay->SetTitle("Ceres", kFALSE);
128
129 cont.Add(fDisplay);
130 }
131
132// const TString oname = Form("reflector%08d.root", fSequence.GetSequence());
133 const TString oname = "ceres.root";
134 return WriteContainer(cont, oname, "RECREATE");
135}
136
137void MJSimulation::SetupHist(MHn &hist) const
138{
139 hist.AddHist("MCorsikaEvtHeader.fTotalEnergy");
140 hist.InitName("Energy");
141 hist.InitTitle("Energy;E [GeV]");
142 hist.SetLog(kTRUE, kTRUE, kFALSE);
143
144 hist.AddHist("MPhotonEvent.GetNumPhotons");
145 hist.InitName("Size");
146 hist.InitTitle("Size;S [#]");
147 hist.SetLog(kTRUE, kTRUE, kFALSE);
148
149 hist.AddHist("MCorsikaEvtHeader.fX/100","MCorsikaEvtHeader.fY/100");
150 hist.SetDrawOption("colz");
151 hist.InitName("Impact;Impact;Impact");
152 hist.InitTitle("Impact;West <--> East [m];South <--> North [m]");
153
154 hist.AddHist("MCorsikaEvtHeader.fFirstInteractionHeight/100000");
155 hist.InitName("Height");
156 hist.InitTitle("FirstInteractionHeight;h [km]");
157
158 hist.AddHist("MCorsikaEvtHeader.fAz*TMath::RadToDeg()", "MCorsikaEvtHeader.fZd*TMath::RadToDeg()");
159 hist.InitName("SkyOrigin;Az;Zd");
160 hist.InitTitle("Sky Origin;Az [deg];Zd [deg]");
161 hist.SetDrawOption("colz");
162
163 TString sin2 = "sin(MCorsikaEvtHeader.fZd)*sin(MCorsikaRunHeader.fZdMin*TMath::DegToRad())";
164 TString cos2 = "cos(MCorsikaEvtHeader.fZd)*cos(MCorsikaRunHeader.fZdMin*TMath::DegToRad())";
165 TString cos = "cos(MCorsikaEvtHeader.fAz-MCorsikaRunHeader.fAzMin*TMath::DegToRad())";
166
167 TString form = "acos("+sin2+"*"+cos+"+"+cos2+")*TMath::RadToDeg()";
168
169 hist.AddHist(form);
170 hist.InitName("ViewCone");
171 hist.InitTitle("Incident Angle;\\alpha [\\deg]");
172}
173
174Bool_t MJSimulation::Process(const MArgs &args)
175{
176 /*
177 if (!fSequence.IsValid())
178 {
179 *fLog << err << "ERROR - Sequence invalid!" << endl;
180 return kFALSE;
181 }
182 */
183
184 //if (!HasWritePermission(GetPathOut()))
185 // return kFALSE;
186
187 *fLog << inf;
188 fLog->Separator(GetDescriptor());
189
190 if (!CheckEnv())
191 return kFALSE;
192
193 *fLog << warn << "FIXME: Monte Carlo simulation: Sequences not supported yet.";
194 //*fLog << fSequence.GetFileName() << endl;
195 *fLog << endl;
196
197 // --------------------------------------------------------------------------------
198
199 //MDirIter iter;
200 //if (fSequence.GetRuns(iter, MSequence::kCalibrated)<=0)
201 // return kFALSE;
202
203 // Setup Parlist
204 MParList plist;
205 plist.AddToList(this); // take care of fDisplay!
206
207 // setup TaskList
208 MTaskList tasks;
209 plist.AddToList(&tasks);
210
211 // --------------------------------------------------------------------------------
212
213 // ---------- FIXME FIXME FIXME -----------
214 // FIXME: Find a better name (What if both cams are identical?)
215 MParEnv env1("GeomAPDs"); // Inheritance!
216 MParEnv env2("MGeomCam"); // Inheritance!
217 env1.SetClassName("MGeomCam");
218 env2.SetClassName("MGeomCam");
219 plist.AddToList(&env1);
220 plist.AddToList(&env2);
221
222 // FIXME: Allow the user to create other reflectors
223 MReflector reflector;
224 plist.AddToList(&reflector);
225
226 plist.FindCreateObj("MPedPhotCam", "MPedPhotFromExtractorRndm");
227
228 MPulseShape shape;
229 plist.AddToList(&shape);
230
231 MRawRunHeader header;
232 header.InitFadcType(3);
233
234 header.SetRunInfo(/*MRawRunHeader::kRTMonteCarlo|*/MRawRunHeader::kRTData, 1, 1);
235 if (args.GetNumArguments()==1)
236 {
237 if (!args.GetArgumentStr(0).CompareTo("pedestal", TString::kIgnoreCase))
238 header.SetRunInfo(/*MRawRunHeader::kRTMonteCarlo|*/MRawRunHeader::kRTPedestal, 1, 1);
239 if (!args.GetArgumentStr(0).CompareTo("calibration", TString::kIgnoreCase))
240 header.SetRunInfo(/*MRawRunHeader::kRTMonteCarlo|*/MRawRunHeader::kRTCalibration, 1, 1);
241 if (!args.GetArgumentStr(0).CompareTo("pointrun", TString::kIgnoreCase))
242 header.SetRunInfo(/*MRawRunHeader::kRTMonteCarlo|*/MRawRunHeader::kRTPointRun, 1, 1);
243 }
244
245 // FIXME: Move to MSimPointingPos, MSimCalibrationSignal
246 header.SetObservation("On", "MonteCarlo");
247 plist.AddToList(&header);
248 // ++++++++ FIXME FIXME FIXME +++++++++++++
249
250 /*
251 MPedestalCam pedcam;
252 pedcam.Init(geomcam.GetNumPixels());
253 for (UInt_t i=0; i<geomcam.GetNumPixels(); i++)
254 pedcam[i].Set(128./header.GetScale(), 22.5/header.GetScale());
255 plist.AddToList(&pedcam);
256 */
257
258 // -------------------------------------------------------------------
259
260 MCorsikaRead read;
261 read.SetForceMode(fForceMode);
262
263 for (int i=0; i<args.GetNumArguments();i ++)
264 read.AddFile(args.GetArgumentStr(i));
265
266 MSimAbsorption absapd("AbsorptionAPDs");
267 MSimAbsorption absmir("AbsorptionMirrors");
268 MSimAbsorption cones("AbsorptionCones");
269 cones.SetUseTheta();
270
271 MSimPointingPos pointing;
272
273 MSimReflector reflect;
274 reflect.SetNameGeomCam("GeomAPDs");
275// MSimStarField stars;
276
277 MContinue cont1("MPhotonEvent.GetNumPhotons<1", "ContEmpty1");
278 MContinue cont2("MPhotonEvent.GetNumPhotons<1", "ContEmpty2");
279
280 // -------------------------------------------------------------------
281
282 MBinning binse( 100, 1, 100000, "BinningEnergy", "log");
283 MBinning binss( 100, 1, 10000000, "BinningSize", "log");
284 MBinning binsi( 50, -250, 250, "BinningImpact");
285 MBinning binsh( 50, 0, 12, "BinningHeight");
286 MBinning binsaz(360, -360, 360, "BinningAz");
287 MBinning binszd( 70, 0, 70, "BinningZd");
288 MBinning binsvc(155, 0, 31, "BinningViewCone");
289 MBinning binstr(150, -25, 125, "BinningTrigger");
290 MBinning binsew(150, 0, 150, "BinningEvtWidth");
291
292 plist.AddToList(&binse);
293 plist.AddToList(&binss);
294 plist.AddToList(&binsi);
295 plist.AddToList(&binsh);
296 plist.AddToList(&binszd);
297 plist.AddToList(&binsaz);
298 plist.AddToList(&binsvc);
299 plist.AddToList(&binstr);
300 plist.AddToList(&binsew);
301
302 MHn mhn1, mhn2, mhn3;
303 SetupHist(mhn1);
304 SetupHist(mhn2);
305 SetupHist(mhn3);
306
307 MH3 mhtp("TriggerPos.fVal-IntendedPulsePos.fVal-MPulseShape.GetPulseWidth");
308 mhtp.SetName("Trigger");
309 mhtp.SetTitle("Trigger position w.r.t. the first photon hitting an APD");
310
311 MH3 mhew("MPhotonStatistics.fTimeLast-MPhotonStatistics.fTimeFirst - MRawRunHeader.GetFreqSampling/1000.*MRawRunHeader.fNumSamplesHiGain - 2*MPulseShape.GetPulseWidth");
312 mhew.SetName("EvtWidth");
313 mhew.SetTitle("Time between first and last photon hitting an APD");
314
315 MFillH fillh1(&mhn1, "", "FillH1");
316 MFillH fillh2(&mhn2, "", "FillH2");
317 MFillH fillh3(&mhn3, "", "FillH3");
318 MFillH filltp(&mhtp, "", "FillTriggerPos");
319 MFillH fillew(&mhew, "", "FillEvtWidth");
320 fillh1.SetNameTab("H1", "Distribution of Muons as simulated");
321 fillh2.SetNameTab("H2", "Distribution of Muons as available after all");
322 fillh3.SetNameTab("H3", "Distribution after trigger");
323 filltp.SetNameTab("TrigPos", "TriggerPosition w.r.t the first photon");
324 fillew.SetNameTab("EvtWidth", "Time between first and last photon hitting an APD");
325
326 MHPhotonEvent planeG(1); // Get from MaxImpact
327 MHPhotonEvent plane0(2); // Get from MReflector
328 MHPhotonEvent plane1(2);
329 MHPhotonEvent plane2(2);
330 MHPhotonEvent plane3(2);
331 MHPhotonEvent plane4(2);
332 MHPhotonEvent planeF1(3); // Get from MGeomCam
333 MHPhotonEvent planeF2(header.IsPointRun()?4:3); // Get from MGeomCam
334
335 //MHPSF psf;
336
337 MFillH fillG(&planeG, "MPhotonEvent", "FillGround");
338 MFillH fill0(&plane0, "MirrorPlane0", "FillReflector");
339 MFillH fill1(&plane1, "MirrorPlane1", "FillCamShadow");
340 MFillH fill2(&plane2, "MirrorPlane2", "FillCandidates");
341 MFillH fill3(&plane3, "MirrorPlane3", "FillReflected");
342 MFillH fill4(&plane4, "MirrorPlane4", "FillFocal");
343 MFillH fillF1(&planeF1, "MPhotonEvent", "FillCamera");
344 MFillH fillF2(&planeF2, "MPhotonEvent", "FillAPDs");
345 //MFillH fillP(&psf, "MPhotonEvent", "FillPSF");
346
347 fillG.SetNameTab("Ground", "Photon distribution at ground");
348 fill0.SetNameTab("Reflector", "Photon distribution at reflector plane w/o camera shadow");
349 fill1.SetNameTab("CamShadow", "Photon distribution at reflector plane w/ camera shadow");
350 fill2.SetNameTab("Candidates", "'Can hit' photon distribution at reflector plane w/ camera shadow");
351 fill3.SetNameTab("Reflected", "Photon distribution after reflector projected to reflector plane");
352 fill4.SetNameTab("Focal", "Photon distribution at focal plane");
353 fillF1.SetNameTab("Camera", "Photon distribution which hit the detector");
354 fillF2.SetNameTab("APDs", "Photon distribution after cones");
355 //fillP.SetNameTab("PSF", "Photon distribution after cones for all mirrors");
356
357 // -------------------------------------------------------------------
358
359 const TString rule1(Form("s/cer([0-9]+)/%s\\/ref$1.root/", Esc(fPathOut).Data()));
360 const TString rule2(Form("s/cer([0-9]+)/%s\\/sig$1.root/", Esc(fPathOut).Data()));
361 const TString rule3(Form("s/cer([0-9]+)/%s\\/cam$1.root/", Esc(fPathOut).Data()));
362
363 MWriteRootFile write1( 2, rule1, fOverwrite?"RECREATE":"NEW", "Reflector file");
364 MWriteRootFile write2( 2, rule2, fOverwrite?"RECREATE":"NEW", "Signal file");
365 MWriteRootFile write3( 2, rule3, fOverwrite?"RECREATE":"NEW", "Camera file");
366
367 write1.SetName("WriteRef");
368 write2.SetName("WriteSig");
369 write3.SetName("WriteCam");
370
371 write1.AddContainer("MCorsikaEvtHeader", "Events", kFALSE);
372 write1.AddContainer("MPhotonEvent", "Events");
373
374 write2.AddContainer("MRawRunHeader", "RunHeaders");
375 write2.AddContainer("MGeomCam", "RunHeaders");
376 write2.AddContainer("MSignalCam", "Events");
377 write2.AddContainer("MRawEvtHeader", "Events");
378 write2.AddContainer("MPedPhotFromExtractorRndm", "RunHeaders");
379 /* Monte Carlo Headers
380 write.AddContainer("MMcRunHeader", "RunHeaders", kFALSE);
381 write.AddContainer("MMcFadcHeader", "RunHeaders", kFALSE);
382 write.AddContainer("MMcTrigHeader", "RunHeaders", kFALSE);
383 write.AddContainer("MMcConfigRunHeader", "RunHeaders", kFALSE);
384 write.AddContainer("MMcCorsikaRunHeader", "RunHeaders", kFALSE);
385 // Monte Carlo
386 write.AddContainer("MMcEvt", "Events", kFALSE);
387 write.AddContainer("MMcTrig", "Events", kFALSE);
388 */// Data tree
389 //write.AddContainer("MPedPhotFromExtractor", "Events");
390 //write.AddContainer("MPedPhotFromExtractorRndm", "Events");
391 //write.AddContainer("MTime", "Events", kFALSE);
392 //write.AddContainer("MRawEvtHeader", "Events");
393 //write.AddContainer("MTriggerPattern", "Events");
394
395 write3.AddContainer("MRawEvtData", "Events");
396 write3.AddContainer("MRawRunHeader", "RunHeaders");
397 write3.AddContainer("MGeomCam", "RunHeaders");
398 write3.AddContainer("MRawEvtHeader", "Events");
399 write3.AddContainer("MPedestalCam", "RunHeaders", kFALSE);
400
401 // -------------------------------------------------------------------
402
403 MGeomApply apply;
404
405 MSimGeomCam simgeom;
406 simgeom.SetNameGeomCam("GeomAPDs");
407
408 MSimCalibrationSignal simcal;
409 simcal.SetNameGeomCam("GeomAPDs");
410
411 // Dark Current: 4MHz = 0.004 GHz
412 // Light of night sky at La Palma: ~ 0.2 ph / cm^2 / sr / ns
413 // 5deg camera: 6e-3 sr
414 // NSB = 0.2*6e-3
415
416 MSimRandomPhotons simnsb;
417// simnsb.SetFreq(0.005, 0.004); // 5MHz/cm^2, 4MHz
418 simnsb.SetFreq(0, 0.04); // 40MHz fixed for each APD
419 simnsb.SetNameGeomCam("GeomAPDs");
420
421 // FIXME: How do we handle star-light? We need the rate but we also
422 // need to process it through the mirror!
423
424 MSimAPD simapd;
425 simapd.SetFreq(0.04); // Must be identical to MSimRanodmPhotons!!
426 simapd.SetNameGeomCam("GeomAPDs");
427
428 MSimExcessNoise simexcnoise;
429 MSimBundlePhotons simsum;
430 MSimSignalCam simsignal;
431 MSimCamera simcam;
432 MSimTrigger simtrig;
433 MSimReadout simdaq;
434
435 MContinue conttrig("TriggerPos.fVal<0", "ContTrigger");
436
437 MParameterD pulpos("IntendedPulsePos");
438// pulpos.SetVal(40); // [ns]
439 plist.AddToList(&pulpos);
440
441 MParameterD trigger("TriggerPos");
442 trigger.SetVal(0);
443 plist.AddToList(&trigger);
444
445 // -------------------------------------------------------------------
446
447 MImgCleanStd clean(7, 4.5);
448 clean.SetMethod(MImgCleanStd::kAbsolute);
449
450 //clean.SetNamePedPhotCam("MPedPhotFromExtractorRndm");
451
452 MHillasCalc hcalc;
453 hcalc.Disable(MHillasCalc::kCalcConc);
454
455 MHCamEvent evt0a(/*10*/0, "Signal", "Average signal;;S [ph]");
456 MHCamEvent evt0d(/*11*/8, "ArrTm", "Time after first photon;;T [ns]");
457 evt0a.SetErrorSpread(kFALSE);
458
459 MFillH fillx0a(&evt0a, "MSignalCam", "FillSignal");
460 MFillH fillx0d(&evt0d, "MSignalCam", "FillArrTm");
461 MFillH fillx1("MHHillas", "MHillas", "FillHillas");
462 //MFillH fillx2("MHHillasExt", "", "FillHillasExt");
463 MFillH fillx3("MHHillasSrc", "MHillasSrc", "FillHillasSrc");
464 //MFillH fillx4("MHImagePar", "MImagePar", "FillImagePar");
465 //MFillH fillx5("MHNewImagePar", "MNewImagePar", "FillNewImagePar");
466
467 // -------------------------------------------------------------------
468
469 if (header.IsDataRun())
470 {
471 tasks.AddToList(&read);
472 tasks.AddToList(&pointing);
473 // if (header.IsPointRun())
474 // tasks.AddToList(&stars);
475 //tasks.AddToList(&print);
476 tasks.AddToList(&fillh1);
477 tasks.AddToList(&fillG);
478 if (!header.IsPointRun())
479 {
480 tasks.AddToList(&absapd);
481 tasks.AddToList(&absmir);
482 }
483 tasks.AddToList(&reflect);
484 if (!header.IsPointRun())
485 {
486 tasks.AddToList(&fill0);
487 tasks.AddToList(&fill1);
488 tasks.AddToList(&fill2);
489 tasks.AddToList(&fill3);
490 tasks.AddToList(&fill4);
491 tasks.AddToList(&fillF1);
492 }
493 tasks.AddToList(&cones);
494 tasks.AddToList(&fillF2);
495 //if (header.IsPointRun())
496 // tasks.AddToList(&fillP);
497 tasks.AddToList(&cont1);
498 if (!header.IsPointRun())
499 tasks.AddToList(&fillh2);
500 }
501 // -------------------------------
502 tasks.AddToList(&apply);
503 if (header.IsDataRun())
504 {
505 tasks.AddToList(&simgeom);
506 tasks.AddToList(&cont2);
507 }
508 if (header.IsPedestalRun() || header.IsCalibrationRun())
509 tasks.AddToList(&simcal);
510 tasks.AddToList(&simnsb);
511 tasks.AddToList(&simapd);
512 tasks.AddToList(&simexcnoise);
513 tasks.AddToList(&simsum);
514 tasks.AddToList(&simcam);
515 if (header.IsDataRun())
516 tasks.AddToList(&simtrig);
517 tasks.AddToList(&conttrig);
518 tasks.AddToList(&simsignal); // What do we do if signal<0?
519 tasks.AddToList(&simdaq);
520 if (!fPathOut.IsNull() && !HasNullOut())
521 {
522 tasks.AddToList(&write1);
523 tasks.AddToList(&write2);
524 tasks.AddToList(&write3);
525 }
526 // -------------------------------
527 if (header.IsDataRun())
528 tasks.AddToList(&fillh3);
529 tasks.AddToList(&filltp);
530 tasks.AddToList(&fillew);
531 tasks.AddToList(&fillx0a);
532// tasks.AddToList(&clean);
533 tasks.AddToList(&hcalc);
534 tasks.AddToList(&fillx0d);
535 tasks.AddToList(&fillx1);
536 //tasks.AddToList(&fillx2);
537 tasks.AddToList(&fillx3);
538 //tasks.AddToList(&fillx4);
539 //tasks.AddToList(&fillx5);
540
541 //-------------------------------------------
542
543 tasks.SetAccelerator(MTask::kAccDontReset|MTask::kAccDontTime);
544
545 // Create and setup the eventloop
546 MEvtLoop evtloop(fName);
547 evtloop.SetParList(&plist);
548 evtloop.SetDisplay(fDisplay);
549 evtloop.SetLogStream(&gLog);
550 if (!SetupEnv(evtloop))
551 return kFALSE;
552
553 // Execute first analysis
554 if (!evtloop.Eventloop(fMaxEvents))
555 {
556 gLog << err << GetDescriptor() << ": Failed." << endl;
557 return kFALSE;
558 }
559
560 //-------------------------------------------
561 // FIXME: Display Spline in tab
562 // FIXME: Display Reflector in tab
563 // FIXME: Display Routing(?) in tab
564 // FIXME: Display Camera(s) in tab
565 //-------------------------------------------
566
567 if (!WriteResult())
568 return kFALSE;
569
570 *fLog << all << GetDescriptor() << ": Done." << endl << endl << endl;;
571
572 return kTRUE;
573}
Note: See TracBrowser for help on using the repository browser.