source: branches/Mars_MC/msimreflector/MSimRays.cc@ 17055

Last change on this file since 17055 was 14153, checked in by tbretz, 12 years ago
Replaced my own homogeneous number rgenerator with gRandom->Circle
File size: 7.4 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:thomas.bretz@epfl.ch>
19!
20! Copyright: CheObs Software Development, 2000-2010
21!
22!
23\* ======================================================================== */
24
25//////////////////////////////////////////////////////////////////////////////
26//
27// MSimRays
28//
29// Task to produce rays from a light source at either infinity or a given
30// height from a given local sky position.
31//
32// The sky position is defined by an MPointingPos object in the parameter
33// list (if none exists, the source is at the reflector axis). Its
34// default name is "MPointingPos".
35//
36// The height of the light/point source is set by SetHeight in units of km.
37// A value <= 0 means infinity.
38//
39// The number of rays produced per event is defined by SetNumPhotons(n).
40// The default is 1000.
41//
42//////////////////////////////////////////////////////////////////////////////
43#include "MSimRays.h"
44
45#include <TMath.h> // root >=5.20
46#include <TRandom.h>
47#include <TRotation.h>
48
49#include "MLog.h"
50#include "MLogManip.h"
51
52#include "MParList.h"
53
54#include "MQuaternion.h"
55
56#include "MPhotonEvent.h"
57#include "MPhotonData.h"
58
59#include "MReflector.h"
60#include "MPointingPos.h"
61
62ClassImp(MSimRays);
63
64using namespace std;
65
66// --------------------------------------------------------------------------
67//
68// Default Constructor.
69//
70MSimRays::MSimRays(const char* name, const char *title)
71 : fEvt(0), fReflector(0), fPointPos(0), fSource(0),
72 fNumPhotons(1000), fHeight(-1),
73 fNameReflector("MReflector"), fNamePointPos("MPointingPos"),
74 fNameSource("Source")
75{
76 fName = name ? name : "MSimRays";
77 fTitle = title ? title : "Task to calculate reflection os a mirror";
78}
79
80// --------------------------------------------------------------------------
81//
82// Search for the necessary parameter containers.
83//
84Int_t MSimRays::PreProcess(MParList *pList)
85{
86 fEvt = (MPhotonEvent*)pList->FindCreateObj("MPhotonEvent");
87 if (!fEvt)
88 return kFALSE;
89
90 if (!pList->FindCreateObj("MCorsikaEvtHeader"))
91 return kFALSE;
92
93 fReflector = (MReflector*)pList->FindObject(fNameReflector, "MReflector");
94 if (!fReflector)
95 {
96 *fLog << inf << fNameReflector << " [MReflector] not found..." << endl;
97 return kFALSE;
98 }
99
100 fSource = (MPointingPos*)pList->FindObject(fNameSource, "MPointingPos");
101 if (!fSource)
102 {
103// *fLog << inf << fNameSource << " [MPointingPos] not found..." << endl;
104// return kFALSE;
105 }
106
107 fPointPos = (MPointingPos*)pList->FindObject(fNamePointPos, "MPointingPos");
108 if (!fPointPos)
109 {
110 *fLog << inf << fNamePointPos << " [MPointingPos] not found..." << endl;
111 return kFALSE;
112 }
113
114 return kTRUE;
115}
116
117// --------------------------------------------------------------------------
118//
119// Converts the photons into the telscope coordinate frame using the
120// pointing position from MPointingPos.
121//
122// Reflects all photons on all mirrors and stores the final photons on
123// the focal plane. Also intermediate photons are stored for debugging.
124//
125Int_t MSimRays::Process()
126{
127 // Get arrays from event container
128 fEvt->Resize(fNumPhotons);
129
130 TClonesArray &arr = fEvt->GetArray();
131
132 const Int_t num = arr.GetEntriesFast();
133
134 const Double_t maxr = fReflector->GetMaxR();
135
136 const Double_t deltazd = fSource ? fSource->GetZdRad() : 0;
137 const Double_t deltaaz = fSource ? fSource->GetAzRad() : 0;
138
139 const Double_t zd = fPointPos->GetZdRad() + deltazd;
140 const Double_t az = fPointPos->GetAzRad() + deltaaz;
141
142
143 // cm -> m
144 // s -> ns
145 // length -> time
146 const Double_t conv = 1./(TMath::C()*100/1e9);
147
148 // Local sky coordinates (direction of telescope axis)
149 //const Double_t zd = fPointing->GetZdRad(); // x==north
150 //const Double_t az = fPointing->GetAzRad();
151
152 // Height of point source [cm] (0 means infinity)
153 const Double_t h = fHeight * 100000;
154
155 // Rotation matrix to derotate sky
156 // For the new coordinate system see the Wiki
157 TRotation rot; // The signs are positive because we align the incident point on ground to the telescope axis
158 rot.RotateX( zd); // Rotate point on ground to align it with the telescope axis
159 rot.RotateZ(-az); // tilt the point from ground to make it parallel to the mirror plane
160
161 Int_t idx = 0;
162 while (idx<num)
163 {
164 MPhotonData &dat = *static_cast<MPhotonData*>(arr.UncheckedAt(idx));
165
166 Double_t x, y;
167 const Double_t r = gRandom->Uniform();
168 gRandom->Circle(x, y, maxr*TMath::Sqrt(r));
169/*
170 Double_t ra = gRandom->Uniform(maxr);
171 Double_t ph = gRandom->Uniform(TMath::TwoPi());
172
173
174 // Get radom incident point on the mirror plane.
175 //const Double_t x = gRandom->Uniform(-maxr, maxr);
176 //const Double_t y = gRandom->Uniform(-maxr, maxr);
177
178 Double_t x = ra*sin(ph);
179 Double_t y = ra*cos(ph);
180
181// if (x*x + y*y > maxr*maxr)
182// continue;
183 */
184 // The is the incident direction of the photon
185 // h==0 means infinitiy
186 const TVector3 u = fHeight>0 ? TVector3(x, y, -h).Unit() : TVector3(0, 0, -1);
187
188 // w is pointing away from the direction the photon comes from
189 // CORSIKA-orig: x(north), y(west), z(up), t(time)
190 // NOW: x(east), y(north), z(up), t(time)
191 MQuaternion p(TVector3(x, y, 0), fHeight>0 ? TMath::Sqrt(x*x + y*y + h*h): 0);
192 MQuaternion w(u, conv);
193
194 // Rotate the coordinates into the reflector's coordinate system.
195 // It is assumed that the z-plane is parallel to the focal plane.
196 // (The reflector coordinate system is defined by the telescope orientation)
197 p *= rot;
198 w *= rot;
199
200 // Now propagate the photon to the z-plane in the new coordinate system
201 p.PropagateZ0(w);
202
203 // Shift the coordinate system to the telescope. Corsika's
204 // coordinate system is always w.r.t. to the particle axis
205 //p += impact;
206
207 // Store new position and direction in the reflector's coordinate frame
208 dat.SetPosition(p);
209 dat.SetDirection(w);
210
211 idx++;
212 }
213
214 // Doesn't seem to be too time consuming. But we could also sort later!
215 // (after cones, inside the camera)
216 // fEvt->Sort(kTRUE);
217
218 return kTRUE;
219}
220
221// --------------------------------------------------------------------------
222//
223// Height: -1
224// NumPhotons: 1000
225//
226Int_t MSimRays::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
227{
228 Bool_t rc = kFALSE;
229 if (IsEnvDefined(env, prefix, "Height", print))
230 {
231 rc = kTRUE;
232 fHeight = GetEnvValue(env, prefix, "Height", fHeight);
233 }
234 if (IsEnvDefined(env, prefix, "NumPhotons", print))
235 {
236 rc = kTRUE;
237 fNumPhotons = GetEnvValue(env, prefix, "NumPhotons", (Int_t)fNumPhotons);
238 }
239
240 return rc;
241}
Note: See TracBrowser for help on using the repository browser.