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

Last change on this file since 18804 was 18451, checked in by ftemme, 10 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.