source: trunk/Mars/mjobs/MJSimReflector.cc@ 19727

Last change on this file since 19727 was 18451, checked in by ftemme, 9 years ago
Added the new program mirrorSimulation which performs the simulation of the reflector. Adapted the Makefiles and added the corresponding mjob MJSimReflector
File size: 13.6 KB
Line 
1/////////////////////////////////////////////////////////////////////////////
2//
3// MJSimReflector
4//
5//
6// This process performs the simulation up to the resulting photons on the camera
7// window:
8//
9// - Simulating of Pointing
10// - Absorption in the atmosphere
11// - Absorption by mirror reflectivity and additional photon acceptance
12// - Simulation of the reflector
13// - Absorption by angular acceptance of winston cones
14// - Absorption by transmission acceptance of winston cones
15// - Absorption by pde of the SiPMs
16//
17// An ASCII output file with the information of all photons is written to disk.
18// The place in the simulation chain can be chosen by changing the fill point of the
19// writeCherPhotonFile into the tasklist
20//
21//
22/////////////////////////////////////////////////////////////////////////////
23
24#include "MJSimReflector.h"
25
26#include <TEnv.h>
27#include <TCanvas.h>
28#include <iostream>
29
30// Core
31#include "MLog.h"
32#include "MLogManip.h"
33
34#include "MArgs.h"
35#include "MDirIter.h"
36#include "MParList.h"
37#include "MTaskList.h"
38#include "MEvtLoop.h"
39
40#include "MStatusDisplay.h"
41
42#include "MWriteAsciiFile.h"
43
44// Tasks
45#include "MCorsikaRead.h"
46#include "MContinue.h"
47#include "MGeomApply.h"
48#include "MParameterCalc.h"
49
50#include "MSimMMCS.h"
51#include "MSimAbsorption.h"
52#include "MSimAtmosphere.h"
53#include "MSimReflector.h"
54#include "MSimPointingPos.h"
55#include "MSimPSF.h"
56#include "MSimGeomCam.h"
57#include "MSimRandomPhotons.h"
58#include "MSimBundlePhotons.h"
59
60// Container
61#include "MRawRunHeader.h"
62#include "MParameters.h"
63#include "MReflector.h"
64#include "MParEnv.h"
65#include "MSpline3.h"
66#include "MParSpline.h"
67#include "MGeomCam.h"
68#include "MMatrix.h"
69
70#include "MPedestalCam.h"
71#include "MPedestalPix.h"
72
73ClassImp(MJSimReflector);
74
75using namespace std;
76
77// --------------------------------------------------------------------------
78//
79// Default constructor.
80//
81// Sets fRuns to 0, fExtractor to NULL, fDataCheck to kFALSE
82//
83MJSimReflector::MJSimReflector(const char *name, const char *title)
84 : fForceMode(kFALSE), fOperationMode(kModeData), fRunNumber(-1)
85{
86 fName = name ? name : "MJSimReflector";
87 fTitle = title ? title : "Standard analysis and reconstruction";
88}
89
90Bool_t MJSimReflector::CheckEnvLocal()
91{
92 fForceMode = GetEnv("ForceMode", fForceMode);
93
94 return kTRUE;
95}
96
97// Setup the header accordingly to the used operation mode
98void MJSimReflector::SetupHeaderOperationMode(MRawRunHeader &header) const
99{
100 switch (fOperationMode)
101 {
102 case kModeData:
103 header.SetRunType(MRawRunHeader::kRTMonteCarlo|MRawRunHeader::kRTData);
104 header.SetRunInfo(0, fRunNumber<0 ? 3 : fRunNumber);
105 break;
106
107 case kModePed:
108 header.SetRunType(MRawRunHeader::kRTMonteCarlo|MRawRunHeader::kRTPedestal);
109 header.SetSourceInfo("Pedestal");
110 header.SetRunInfo(0, fRunNumber<0 ? 1 : fRunNumber);
111 break;
112
113 case kModeCal:
114 header.SetRunType(MRawRunHeader::kRTMonteCarlo|MRawRunHeader::kRTCalibration);
115 header.SetSourceInfo("Calibration");
116 header.SetRunInfo(0, fRunNumber<0 ? 2 : fRunNumber);
117 break;
118
119 case kModePointRun:
120 header.SetRunType(MRawRunHeader::kRTMonteCarlo|MRawRunHeader::kRTPointRun);
121 header.SetRunInfo(0, fRunNumber<0 ? 0 : fRunNumber);
122 break;
123 }
124}
125
126
127Bool_t MJSimReflector::Process(const MArgs &args, const MSequence &seq)
128{
129 // --------------------------------------------------------------------------------
130 // --------------------------------------------------------------------------------
131 // Initialization
132 // --------------------------------------------------------------------------------
133 // --------------------------------------------------------------------------------
134
135 // --------------------------------------------------------------------------------
136 // Logging output:
137 // --------------------------------------------------------------------------------
138 // - description of the job itself
139 // - list of the
140 *fLog << inf;
141 fLog->Separator(GetDescriptor());
142
143 if (!CheckEnv())
144 return kFALSE;
145
146 if (seq.IsValid())
147 *fLog << fSequence.GetFileName() << endl;
148 else
149 *fLog << "Input: " << args.GetNumArguments() << "-files" << endl;
150 *fLog << endl;
151
152 MDirIter iter;
153 if (seq.IsValid() && seq.GetRuns(iter, MSequence::kCorsika)<=0)
154 {
155 *fLog << err << "ERROR - Sequence valid but without files." << endl;
156 return kFALSE;
157 }
158 // --------------------------------------------------------------------------------
159 // Setup MParList and MTasklist
160 // --------------------------------------------------------------------------------
161 MParList plist;
162 plist.AddToList(this); // take care of fDisplay!
163 // setup TaskList
164 MTaskList tasks;
165 plist.AddToList(&tasks);
166
167 // --------------------------------------------------------------------------------
168 // --------------------------------------------------------------------------------
169 // Parameter Container Setup
170 // --------------------------------------------------------------------------------
171 // --------------------------------------------------------------------------------
172
173 // --------------------------------------------------------------------------------
174 // Setup container for the reflector, the cones and the camera geometry
175 // --------------------------------------------------------------------------------
176 MParEnv env0("Reflector");
177 MParEnv env1("GeomCones");
178 MParEnv env2("MGeomCam");
179 env0.SetClassName("MReflector");
180 env1.SetClassName("MGeomCam");
181 env2.SetClassName("MGeomCam");
182 plist.AddToList(&env0);
183 plist.AddToList(&env1);
184 plist.AddToList(&env2);
185 // --------------------------------------------------------------------------------
186 // Setup spline containers for:
187 // --------------------------------------------------------------------------------
188 // - the different photon acceptance splines
189 MParSpline splinepde("PhotonDetectionEfficiency");
190 MParSpline splinemirror("MirrorReflectivity");
191 MParSpline splinecones("ConesAngularAcceptance");
192 MParSpline splinecones2("ConesTransmission");
193 MParSpline splineAdditionalPhotonAcceptance("AdditionalPhotonAcceptance");
194 plist.AddToList(&splinepde);
195 plist.AddToList(&splinemirror);
196 plist.AddToList(&splinecones);
197 plist.AddToList(&splinecones2);
198 plist.AddToList(&splineAdditionalPhotonAcceptance);
199
200 // --------------------------------------------------------------------------------
201 // Setup container for the MC Run Header and the Raw Run Header
202 // --------------------------------------------------------------------------------
203 plist.FindCreateObj("MMcRunHeader");
204
205 MRawRunHeader header;
206 header.SetValidMagicNumber();
207 SetupHeaderOperationMode(header);
208 header.SetObservation("On", "MonteCarlo");
209 plist.AddToList(&header);
210
211 // --------------------------------------------------------------------------------
212 // --------------------------------------------------------------------------------
213 // Tasks Setup
214 // --------------------------------------------------------------------------------
215 // --------------------------------------------------------------------------------
216
217 // --------------------------------------------------------------------------------
218 // Corsika read setup
219 // --------------------------------------------------------------------------------
220 MCorsikaRead read;
221 read.SetForceMode(fForceMode);
222
223 if (!seq.IsValid())
224 {
225 for (int i=0; i<args.GetNumArguments(); i++)
226 read.AddFile(args.GetArgumentStr(i));
227 }
228 else
229 read.AddFiles(iter);
230
231 // --------------------------------------------------------------------------------
232 // Precut after reading
233 // --------------------------------------------------------------------------------
234 // I am not sure if here there is a default for numPhotons < 10
235 MContinue precut("", "PreCut");
236 precut.IsInverted();
237 precut.SetAllowEmpty();
238
239 // --------------------------------------------------------------------------------
240 // Simmmcs and calculation of the incident angle.
241 // --------------------------------------------------------------------------------
242 MSimMMCS simmmcs;
243
244 // The different strings defines the calculation of the incident angle
245 const TString sin2 = "sin(MCorsikaEvtHeader.fZd)*sin(MCorsikaRunHeader.fZdMin*TMath::DegToRad())";
246 const TString cos2 = "cos(MCorsikaEvtHeader.fZd)*cos(MCorsikaRunHeader.fZdMin*TMath::DegToRad())";
247 const TString cos = "cos(MCorsikaEvtHeader.fAz-MCorsikaRunHeader.fAzMin*TMath::DegToRad())";
248 const TString form = "acos("+sin2+"*"+cos+"+"+cos2+")*TMath::RadToDeg()";
249
250 MParameterCalc calcangle(form, "CalcIncidentAngle");
251 calcangle.SetNameParameter("IncidentAngle");
252
253 // --------------------------------------------------------------------------------
254 // Setup of the different Absorptions
255 // --------------------------------------------------------------------------------
256 // - atmosphere
257 // - PDE of the SiPMs
258 // - mirror reflectivity
259 // - cones angular acceptance
260 // - cones transmission
261 // - additional photon acceptance
262 MSimAtmosphere simatm;
263 MSimAbsorption absapd("SimPhotonDetectionEfficiency");
264 MSimAbsorption absmir("SimMirrorReflectivity");
265 MSimAbsorption cones("SimConesAngularAcceptance");
266 MSimAbsorption cones2("SimConesTransmission");
267 MSimAbsorption additionalPhotonAcceptance("SimAdditionalPhotonAcceptance");
268 absapd.SetParName("PhotonDetectionEfficiency");
269 absmir.SetParName("MirrorReflectivity");
270 cones.SetParName("ConesAngularAcceptance");
271 cones.SetUseTheta();
272 cones2.SetParName("ConesTransmission");
273 additionalPhotonAcceptance.SetParName("AdditionalPhotonAcceptance");
274 // --------------------------------------------------------------------------------
275 // Pointing position and reflector simulation
276 // --------------------------------------------------------------------------------
277 MSimPointingPos pointing;
278 MSimReflector reflect;
279 reflect.SetNameGeomCam("GeomCones");
280 reflect.SetNameReflector("Reflector");
281 // --------------------------------------------------------------------------------
282 // Seupt of some continue conditions
283 // --------------------------------------------------------------------------------
284 MContinue cont1("MPhotonEvent.GetNumPhotons<1", "ContEmpty1");
285 MContinue cont2("MPhotonEvent.GetNumPhotons<1", "ContEmpty2");
286 MContinue cont3("MPhotonEvent.GetNumPhotons<2", "ContEmpty3");
287 cont1.SetAllowEmpty();
288 cont2.SetAllowEmpty();
289 cont3.SetAllowEmpty();
290
291 // --------------------------------------------------------------------------------
292 // Geom apply, Simulation psf, and simulation of the pixel geometry
293 // --------------------------------------------------------------------------------
294 // Be awere MGeomApply only resizes parameter container which heritates from
295 // MGeomCam to the actual size of the camera
296 MGeomApply apply;
297
298 // individual single mirror psf
299 MSimPSF simpsf;
300
301 // Simulate the geometry of the pixels (so which photon hit which pixel)
302 MSimGeomCam simgeom;
303 simgeom.SetNameGeomCam("GeomCones");
304
305 // --------------------------------------------------------------------------------
306 // Setup of the Write Task
307 // --------------------------------------------------------------------------------
308 // Setup of the outputpath
309 if (!fFileOut.IsNull())
310 {
311 const Ssiz_t dot = fFileOut.Last('.');
312 const Ssiz_t slash = fFileOut.Last('/');
313 if (dot>slash)
314 fFileOut = fFileOut.Remove(dot);
315 }
316
317 const char *fmt = Form("%s/%08d%%s.csv", fPathOut.Data(), header.GetRunNumber());
318 TString outputFilePath(Form(fmt, Form("%s", "" )));
319
320 MWriteAsciiFile writePhotonFile(outputFilePath,"MPhotonEvent");
321
322 // --------------------------------------------------------------------------------
323 // --------------------------------------------------------------------------------
324 // Filling of tasks in tasklist
325 // --------------------------------------------------------------------------------
326 // --------------------------------------------------------------------------------
327
328 if (header.IsDataRun())
329 {
330 tasks.AddToList(&read);
331 tasks.AddToList(&pointing);
332 tasks.AddToList(&simmmcs);
333 tasks.AddToList(&simatm);
334 tasks.AddToList(&calcangle);
335 if (!header.IsPointRun())
336 {
337 tasks.AddToList(&absmir);
338 tasks.AddToList(&additionalPhotonAcceptance);
339 }
340 tasks.AddToList(&reflect);
341 }
342 tasks.AddToList(&apply);
343 if (header.IsDataRun())
344 {
345// tasks.AddToList(&simpsf);
346 tasks.AddToList(&simgeom);
347 tasks.AddToList(&writePhotonFile);
348 tasks.AddToList(&cones);
349 tasks.AddToList(&cones2);
350 if (!header.IsPointRun())
351 {
352 tasks.AddToList(&absapd);
353 }
354 }
355 tasks.SetAccelerator(MTask::kAccDontReset|MTask::kAccDontTime);
356
357 // Create and setup the eventloop
358 MEvtLoop evtloop(fName);
359 evtloop.SetParList(&plist);
360 evtloop.SetDisplay(fDisplay);
361 evtloop.SetLogStream(&gLog);
362 if (!SetupEnv(evtloop))
363 return kFALSE;
364
365 header.Print();
366
367 // Execute first analysis
368 if (!evtloop.Eventloop(fMaxEvents))
369 {
370 gLog << err << GetDescriptor() << ": Failed." << endl;
371 return kFALSE;
372 }
373
374// if (!WriteResult(plist, seq, header.GetRunNumber()))
375// return kFALSE;
376
377 *fLog << all << GetDescriptor() << ": Done." << endl << endl << endl;;
378
379 return kTRUE;
380}
Note: See TracBrowser for help on using the repository browser.