source: trunk/MagicSoft/Mars/msimreflector/MSimReflector.cc@ 9456

Last change on this file since 9456 was 9356, checked in by tbretz, 16 years ago
*** empty log message ***
File size: 20.3 KB
Line 
1/* ======================================================================== *\
2!
3! *
4! * This file is part of CheObs, the Modular 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 appears 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: CheObs Software Development, 2000-2009
21!
22!
23\* ======================================================================== */
24
25//////////////////////////////////////////////////////////////////////////////
26//
27// MSimReflector
28//
29// fDetectorMargin is a margin (in mm) which is given to the
30// MGeomCam::HitDetector. It should define a margin around the area
31// defined in HitDetector on the focal plane in which photons are kept.
32// Usually this can be 0 because photons not hitting the detector are
33// obsolete except they can later be "moved" inside the detector, e.g.
34// if you use MSimPSF to emulate a PSF by moving photons randomly
35// on the focal plane.
36//
37//////////////////////////////////////////////////////////////////////////////
38#include "MSimReflector.h"
39
40#include <TMath.h>
41#include <TRandom.h>
42
43#include "MGeomCam.h"
44
45#include "MLog.h"
46#include "MLogManip.h"
47
48#include "MParList.h"
49
50#include "MQuaternion.h"
51#include "MMirror.h"
52#include "MReflector.h"
53#include "MReflection.h"
54
55#include "MCorsikaEvtHeader.h"
56//#include "MCorsikaRunHeader.h"
57
58#include "MPhotonEvent.h"
59#include "MPhotonData.h"
60
61#include "MPointingPos.h"
62
63ClassImp(MSimReflector);
64
65using namespace std;
66
67// USEFUL CORSIKA OPTIONS:
68// NOCLONG
69
70// --------------------------------------------------------------------------
71//
72// Default Constructor.
73//
74MSimReflector::MSimReflector(const char* name, const char *title)
75 : fEvt(0), fMirror0(0), fMirror1(0), fMirror2(0), fMirror3(0),
76 fMirror4(0), /*fRunHeader(0),*/ fEvtHeader(0), fReflector(0),
77 fGeomCam(0), fPointing(0), fDetectorMargin(0)
78{
79 fName = name ? name : "MSimReflector";
80 fTitle = title ? title : "Task to calculate reflection os a mirror";
81}
82
83// --------------------------------------------------------------------------
84//
85// Search for the necessary parameter containers.
86//
87Int_t MSimReflector::PreProcess(MParList *pList)
88{
89 fMirror0 = (MPhotonEvent*)pList->FindCreateObj("MPhotonEvent", "MirrorPlane0");
90 if (!fMirror0)
91 return kFALSE;
92 fMirror1 = (MPhotonEvent*)pList->FindCreateObj("MPhotonEvent", "MirrorPlane1");
93 if (!fMirror1)
94 return kFALSE;
95 fMirror2 = (MPhotonEvent*)pList->FindCreateObj("MPhotonEvent", "MirrorPlane2");
96 if (!fMirror2)
97 return kFALSE;
98 fMirror3 = (MPhotonEvent*)pList->FindCreateObj("MPhotonEvent", "MirrorPlane3");
99 if (!fMirror3)
100 return kFALSE;
101 fMirror4 = (MPhotonEvent*)pList->FindCreateObj("MPhotonEvent", "MirrorPlane4");
102 if (!fMirror4)
103 return kFALSE;
104
105 fReflector = (MReflector*)pList->FindObject("Reflector", "MReflector");
106 if (!fReflector)
107 {
108 *fLog << err << "Reflector [MReflector] not found... aborting." << endl;
109 return kFALSE;
110 }
111
112 if (fReflector->GetNumMirrors()==0)
113 {
114 *fLog << err << "ERROR - Reflector doesn't contain a single mirror." << endl;
115 return kFALSE;
116 }
117
118 fGeomCam = (MGeomCam*)pList->FindObject(fNameGeomCam, "MGeomCam");
119 if (!fGeomCam)
120 {
121 *fLog << inf << fNameGeomCam << " [MGeomCam] not found..." << endl;
122
123 fGeomCam = (MGeomCam*)pList->FindObject("MGeomCam");
124 if (!fGeomCam)
125 {
126 *fLog << err << "MGeomCam not found... aborting." << endl;
127 return kFALSE;
128 }
129 }
130
131 fEvt = (MPhotonEvent*)pList->FindObject("MPhotonEvent");
132 if (!fEvt)
133 {
134 *fLog << err << "MPhotonEvent not found... aborting." << endl;
135 return kFALSE;
136 }
137 /*
138 fRunHeader = (MCorsikaRunHeader*)pList->FindObject("MCorsikaRunHeader");
139 if (!fRunHeader)
140 {
141 *fLog << err << "MCorsikaRunHeader not found... aborting." << endl;
142 return kFALSE;
143 }
144 */
145 fEvtHeader = (MCorsikaEvtHeader*)pList->FindObject("MCorsikaEvtHeader");
146 if (!fEvtHeader)
147 {
148 *fLog << err << "MCorsikaEvtHeader not found... aborting." << endl;
149 return kFALSE;
150 }
151
152 fPointing = (MPointingPos*)pList->FindObject(/*"PointingCorsika",*/ "MPointingPos");
153 if (!fPointing)
154 {
155 *fLog << err << "MPointingPos not found... aborting." << endl;
156 return kFALSE;
157 }
158
159 return kTRUE;
160}
161
162// --------------------------------------------------------------------------
163//
164// The main point of calculating the reflection is to determine the
165// coincidence point of the particle trajectory on the mirror surface.
166//
167// If the position and the trajectory of a particle is known it is enough
168// to calculate the z-value of coincidence. x and y are then well defined.
169//
170// Since the problem (mirror) has a rotational symmetry we only have to care
171// about the distance from the z-axis.
172//
173// Given:
174//
175// p: position vector of particle (z=0)
176// u: direction vector of particle
177// F: Focal distance of the mirror
178//
179// We define:
180//
181// q := (px, py )
182// v := (ux/uz, uy/uz)
183// r^2 := x^2 + y^2
184//
185//
186// Distance from z-axis:
187// ---------------------
188//
189// q' = q - z*v (z>0)
190//
191// Calculate distance r (|q|)
192//
193// r^2 = (px-z*ux)^2 + (py-z*uy)^2
194// r^2 = px^2+py^2 + z^2*(ux^2+uy^2) - 2*z*(px*ux+py*uy)
195// r^2 = |q|^2 + z^2*|v|^2 - 2*z* q*v
196//
197//
198// Spherical Mirror Surface: (distance of surface point from 0/0/0)
199// -------------------------
200//
201// Sphere: r^2 + z^2 = R^2 | Parabola: z = p*r^2
202// Mirror: r^2 + (z-R)^2 = R^2 | Mirror: z = p*r^2
203// |
204// Focal length: F=R/2 | Focal length: F = 1/4p
205// |
206// r^2 + (z-2*F)^2 = (2*F)^2 | z = F/4*r^2
207// |
208// z = -sqrt(4*F*F - r*r) + 2*F |
209// z-2*F = -sqrt(4*F*F - r*r) |
210// (z-2*F)^2 = 4*F*F - r*r |
211// z^2-4*F*z+4*F^2 = 4*F*F - r*r (4F^2-r^2>0) | z - F/4*r^2 = 0
212// z^2-4*F*z+r^2 = 0
213//
214// Find the z for which our particle has the same distance from the z-axis
215// as the mirror surface.
216//
217// substitute r^2
218//
219//
220// Equation to solve:
221// ------------------
222//
223// z^2*(1+|v|^2) - 2*z*(2*F+q*v) + |q|^2 = 0 | z^2*|v|^2 - 2*(2/F+q*v)*z + |q|^2 = 0
224//
225// z = (-b +- sqrt(b*b - 4ac))/(2*a)
226//
227// a = 1+|v|^2 | a = |v|^2
228// b = - 2*(2*F+q*v) | b = - 2*(2/F+q*v)
229// c = |q|^2 | c = |q|^2
230// |
231//
232// substitute b := 2*b'
233//
234// z = (-2*b' +- 2*sqrt(b'*b' - ac))/(2*a)
235// z = (- b' +- sqrt(b'*b' - ac))/a
236// z = (-b'/a +- sqrt(b'*b' - ac))/a
237//
238// substitute f := b'/a
239//
240// z = (-f +- sqrt(f^2 - c/a)
241//
242// =======================================================================================
243//
244// After z of the incident point has been determined the position p is
245// propagated along u to the plane with z=z. Now it is checked if the
246// mirror was really hit (this is implemented in HasHit).
247// From the position on the surface and the mirrors curvature we can
248// now calculate the normal vector at the incident point.
249// This normal vector is smeared out with MMirror::PSF (basically a
250// random gaussian) and then the trajectory is reflected on the
251// resulting normal vector.
252//
253Bool_t MMirror::ExecuteReflection(MQuaternion &p, MQuaternion &u) const
254{
255 // If the z-componenet of the direction vector is normalized to 1
256 // the calculation of the incident points becomes very simple and
257 // the resulting z is just the z-coordinate of the incident point.
258 const TVector2 v(u.XYvector()/u.Z());
259 const TVector2 q(p.XYvector());
260
261 // Radius of curvature
262 const Double_t G = 2*fFocalLength;
263
264 // Find the incident point of the vector to the mirror
265 // u corresponds to downwaqrd going particles, thus we use -u here
266 const Double_t b = G - q*v;
267 const Double_t a = v.Mod2();
268 const Double_t c = q.Mod2();
269
270 // Solution for q spherical (a+1) (parabolic mirror (a) instead of (a+1))
271 const Double_t f = b/(a+1);
272 const Double_t g = c/(a+1);
273
274 // Solution of second order polynomial (transformed: a>0)
275 // (The second solution can be omitted, it is the intersection
276 // with the upper part of the sphere)
277 // const Double_t dz = a==0 ? c/(2*b) : f - TMath::Sqrt(f*f - g);
278 const Double_t z = f - TMath::Sqrt(f*f - g);
279
280 // Move the photon along its trajectory to the x/y plane of the
281 // mirror's coordinate frame. Therefor stretch the vector
282 // until its z-component is the distance from the vector origin
283 // until the vector hits the mirror surface.
284 // p += z/u.Z()*u;
285 // p is at the mirror plane and we want to propagate back to the mirror surface
286 p.PropagateZ(u, z);
287
288 // MirrorShape: Now check if the photon really hit the mirror or just missed it
289 if (!HasHit(p))
290 return kFALSE;
291
292 // Get normal vector for reflection by calculating the derivatives
293 // of a spherical mirror along x and y
294 const Double_t d = TMath::Sqrt(G*G - p.R2());
295
296 // This is a normal vector at the incident point
297 TVector3 n(p.X(), p.Y(), -d);
298 // This is the obvious solution for the normal vector
299 // TVector3 n(-p.X()/d, -p.Y()/d, 1));
300
301 if (fSigmaPSF>0)
302 n += SimPSF(n, fFocalLength, fSigmaPSF/10); // Convert from mm to cm
303
304 // Changes also the sign of the z-direction of flight
305 // This is faster giving identical results
306 u *= MReflection(n);
307 //u *= MReflection(p.X(), p.Y(), -d);
308
309 return kTRUE;
310}
311
312// --------------------------------------------------------------------------
313//
314// Converts the coordinates into the coordinate frame of the mirror.
315// Executes the reflection calling ExecuteReflection and converts
316// the coordinates back.
317// Depending on whether the mirror was hit kTRUE or kFALSE is returned.
318// It the mirror was not hit the result coordinates are wrong.
319//
320Bool_t MMirror::ExecuteMirror(MQuaternion &p, MQuaternion &u) const
321{
322 // Move the mirror to the point of origin and rotate the position into
323 // the individual mirrors coordinate frame.
324 // Rotate the direction vector into the mirror's coordinate frame
325 p -= fPos;
326 p *= fTilt;
327 u *= fTilt;
328
329 // Move the photon along its trajectory to the x/y plane of the
330 // mirror's coordinate frame. Therefor stretch the vector
331 // until its z-component vanishes.
332 //p -= p.Z()/u.Z()*u;
333
334 // p is at the reflector plane and we want to propagate back to the mirror plane
335 p.PropagateZ0(u);
336
337 // Now try to propagate the photon from the plane to the mirror
338 // and reflect its direction vector on the mirror.
339 if (!ExecuteReflection(p, u))
340 return kFALSE;
341
342 // Derotate from mirror coordinates and shift the photon back to
343 // reflector coordinates.
344 // Derotate the direction vector
345 u *= fTilt.Inverse();
346 p *= fTilt.Inverse();
347 p += fPos;
348
349 return kTRUE;
350}
351
352// Jeder Spiegel sollte eine Liste aller andern Spiegel in der
353// reihenfolge Ihrer Entfernung enthalten. Wir starten mit der Suche
354// immer beim zuletzt getroffenen Spiegel!
355//
356// --------------------------------------------------------------------------
357//
358// Loops over all mirrors of the reflector. After doing a rough check
359// whether the mirror can be hit at all the reflection is executed
360// calling the ExecuteMirror function of the mirrors.
361//
362// If a mirror was hit its index is retuened, -1 otherwise.
363//
364// FIXME: Do to lopping over all mirrors for all photons this is the
365// most time consuming function in teh reflector simulation. By a more
366// intelligent way of finding the right mirror then just testing all
367// this could be accelerated a lot.
368//
369Int_t MReflector::ExecuteReflector(MQuaternion &p, MQuaternion &u) const
370{
371 //static const TObjArray *arr = &((MMirror*)fMirrors[0])->fNeighbors;
372
373 // This way of access is somuch faster than the program is
374 // a few percent slower if accessed by UncheckedAt
375 const MMirror **s = GetFirstPtr();
376 const MMirror **e = s+GetNumMirrors();
377 //const MMirror **s = (const MMirror**)fMirrors.GetObjectRef(0);
378 //const MMirror **e = s+fMirrors.GetEntriesFast();
379 //const MMirror **s = (const MMirror**)arr->GetObjectRef(0);
380 //const MMirror **e = s+arr->GetEntriesFast();
381
382 // Loop over all mirrors
383 for (const MMirror **m=s; m<e; m++)
384 {
385 const MMirror &mirror = **m;
386
387 // FIXME: Can we speed up using lookup tables or
388 // indexed tables?
389
390 // MirrorShape: Check if this mirror can be hit at all
391 // This is to avoid time consuming calculation if there is no
392 // chance of a coincidence.
393 // FIXME: Inmprove search algorithm (2D Binary search?)
394 if (!mirror.CanHit(p))
395 continue;
396
397 // Make a local copy of position and direction which can be
398 // changed by ExecuteMirror.
399 MQuaternion q(p);
400 MQuaternion v(u);
401
402 // Check if this mirror is hit, and if it is hit return
403 // the reflected position and direction vector.
404 // If the mirror is missed we go on with the next mirror.
405 if (!mirror.ExecuteMirror(q, v))
406 continue;
407
408 // We hit a mirror. Restore the local copy of position and
409 // direction back into p und u.
410 p = q;
411 u = v;
412
413 //arr = &mirror->fNeighbors;
414
415 return m-s;
416 }
417
418 return -1;
419}
420
421// --------------------------------------------------------------------------
422//
423// Converts the photons into the telscope coordinate frame using the
424// pointing position from MPointingPos.
425//
426// Reflects all photons on all mirrors and stores the final photons on
427// the focal plane. Also intermediate photons are stored for debugging.
428//
429Int_t MSimReflector::Process()
430{
431 // Get arrays from event container
432 TClonesArray &arr = fEvt->GetArray();
433
434 // Because we knwo in advance what the maximum storage space could
435 // be we allocated it in advance (or shrink it if it was extremely
436 // huge before)
437 // Note, that the drawback is that an extremly large event
438 // will take about five times its storage space
439 // for a moment even if a lot from it is unused.
440 // It will be freed in the next step.
441 fMirror0->Resize(arr.GetEntriesFast()); // Free memory of allocated MPhotonData
442 fMirror2->Resize(arr.GetEntriesFast()); // Free memory of allocated MPhotonData
443 fMirror3->Resize(arr.GetEntriesFast()); // Free memory of allocated MPhotonData
444 fMirror4->Resize(arr.GetEntriesFast()); // Free memory of allocated MPhotonData
445
446 // Initialize mirror properties
447 const Double_t F = fGeomCam->GetCameraDist()*100; // Focal length [cm]
448
449 // Local sky coordinates (direction of telescope axis)
450 const Double_t zd = fPointing->GetZdRad(); // x==north
451 const Double_t az = fPointing->GetAzRad();
452
453 // Rotation matrix to derotate sky
454 // For the new coordinate system see the Wiki
455 TRotation rot; // The signs are positive because we align the incident point on ground to the telescope axis
456 rot.RotateZ( az); // Rotate point on ground to align it with the telescope axis
457 rot.RotateX(-zd); // tilt the point from ground to make it parallel to the mirror plane
458
459 // Now get the impact point from Corsikas output
460 const TVector3 impact(fEvtHeader->GetX(), fEvtHeader->GetY(), 0);
461
462 // Counter for number of total and final events
463 UInt_t cnt[6] = { 0, 0, 0, 0, 0, 0 };
464
465 const Int_t num = arr.GetEntriesFast();
466 for (Int_t idx=0; idx<num; idx++)
467 {
468 MPhotonData *dat = static_cast<MPhotonData*>(arr.UncheckedAt(idx));
469
470 // w is pointing away from the direction the photon comes from
471 // CORSIKA-orig: x(north), y(west), z(up), t(time)
472 // NOW: x(east), y(north), z(up), t(time)
473 MQuaternion p(dat->GetPosQ()); // z=0
474 MQuaternion w(dat->GetDirQ()); // z<0
475
476 // Shift the coordinate system to the telescope. Corsika's
477 // coordinate system is always w.r.t. to the particle axis
478 p -= impact;
479
480 // Rotate the coordinates into the reflector's coordinate system.
481 // It is assumed that the z-plane is parallel to the focal plane.
482 // (The reflector coordinate system is defined by the telescope orientation)
483 p *= rot;
484 w *= rot;
485
486 // ---> Simulate star-light!
487 // w.fVectorPart.SetXYZ(0.2/17, 0.2/17, -(1-TMath::Hypot(0.3, 0.2)/17));
488
489 // Now propagate the photon to the z-plane in the new coordinate system
490 p.PropagateZ0(w);
491
492 // Store new position and direction in the reflector's coordinate frame
493 dat->SetPosition(p);
494 dat->SetDirection(w);
495
496 (*fMirror0)[cnt[0]++] = *dat;
497 //*static_cast<MPhotonData*>(cpy0.UncheckedAt(cnt[0]++)) = *dat;
498
499 // Check if the photon has hit the camera housing and holding
500 if (fGeomCam->HitFrame(p, w))
501 continue;
502
503 // FIXME: Do we really need this one??
504 //(*fMirror1)[cnt[1]++] = *dat;
505 //*static_cast<MPhotonData*>(cpy1.UncheckedAt(cnt[1]++)) = *dat;
506
507 // Check if the reflector can be hit at all
508 if (!fReflector->CanHit(p))
509 continue;
510
511 (*fMirror2)[cnt[2]++] = *dat;
512 //*static_cast<MPhotonData*>(cpy2.UncheckedAt(cnt[2]++)) = *dat;
513
514 // Now execute the reflection of the photon on the mirrors' surfaces
515 const Int_t num = fReflector->ExecuteReflector(p, w);
516 if (num<0)
517 continue;
518
519 // Set new position and direction (w.r.t. to the reflector's coordinate system)
520 // Set also the index of the mirror which was hit as tag.
521 dat->SetTag(num);
522 dat->SetPosition(p);
523 dat->SetDirection(w);
524
525 (*fMirror3)[cnt[3]++] = *dat;
526 //*static_cast<MPhotonData*>(cpy3.UncheckedAt(cnt[3]++)) = *dat;
527
528 // Propagate the photon along its trajectory to the focal plane z=F
529 p.PropagateZ(w, F);
530
531 // Store new position
532 dat->SetPosition(p);
533
534 (*fMirror4)[cnt[4]++] = *dat;
535 //*static_cast<MPhotonData*>(cpy4.UncheckedAt(cnt[4]++)) = *dat;
536
537 // FIXME: It make make sense to move this out of this class
538 // It is detector specific not reflector specific
539 // Discard all photons which definitly can not hit the detector surface
540 if (!fGeomCam->HitDetector(p, fDetectorMargin))
541 continue;
542
543 // Copy this event to the next 'new' in the list
544 *static_cast<MPhotonData*>(arr.UncheckedAt(cnt[5]++)) = *dat;
545 }
546
547 // Now we shrink the array to a storable size (for details see
548 // MPhotonEvent::Shrink).
549 fMirror0->Shrink(cnt[0]);
550 //fMirror1->Shrink(cnt[1]);
551 fMirror2->Shrink(cnt[2]);
552 fMirror3->Shrink(cnt[3]);
553 fMirror4->Shrink(cnt[4]);
554 fEvt->Shrink(cnt[5]);
555
556 // Doesn't seem to be too time consuming. But we could also sort later!
557 // (after cones, inside the camera)
558 fEvt->Sort(kTRUE);
559
560 // FIXME FIXME FIXME: Set maxindex, first and last time.
561 // SetMaxIndex(fReflector->GetNumMirrors()-1)
562 // if (fEvt->GetNumPhotons())
563 // {
564 // SetTime(fEvt->GetFirst()->GetTime(), fEvt->GetLast()->GetTime());
565 // }
566
567 return kTRUE;
568}
569
570// --------------------------------------------------------------------------
571//
572// DetectorMargin: 0
573//
574Int_t MSimReflector::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
575{
576 Bool_t rc = kFALSE;
577 if (IsEnvDefined(env, prefix, "DetectorMargin", print))
578 {
579 rc = kTRUE;
580 fDetectorMargin = GetEnvValue(env, prefix, "DetectorMargin", 0);
581 }
582
583 return rc;
584}
Note: See TracBrowser for help on using the repository browser.