source: tags/Mars-V0.9.1/mimage/MStereoPar.cc

Last change on this file was 2770, checked in by moralejo, 21 years ago
*** empty log message ***
File size: 10.8 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): Abelardo Moralejo 11/2003 <mailto:moralejo@pd.infn.it>
19!
20! Copyright: MAGIC Software Development, 2000-2003
21!
22!
23\* ======================================================================== */
24
25/////////////////////////////////////////////////////////////////////////////
26//
27// MStereoPar
28//
29// Storage Container for shower parameters estimated using the information
30// from two telescopes (presently for MC studies)
31//
32//
33/////////////////////////////////////////////////////////////////////////////
34#include "MStereoPar.h"
35#include <fstream>
36
37#include "MLog.h"
38#include "MLogManip.h"
39
40#include "MHillas.h"
41#include "MMcEvt.hxx"
42#include "MGeomCam.h"
43
44
45ClassImp(MStereoPar);
46
47using namespace std;
48
49// --------------------------------------------------------------------------
50//
51// Default constructor.
52//
53MStereoPar::MStereoPar(const char *name, const char *title)
54{
55 fName = name ? name : "MStereoPar";
56 fTitle = title ? title : "Stereo image parameters";
57
58
59}
60
61
62// --------------------------------------------------------------------------
63//
64void MStereoPar::Reset()
65{
66 fCoreX = 0.;
67 fCoreY = 0.;
68 fSourceX = 0.;
69 fSourceY = 0.;
70}
71
72
73// --------------------------------------------------------------------------
74//
75// Calculation of shower parameters
76//
77void MStereoPar::Calc(const MHillas &hillas1, const MMcEvt &mcevt1, const MGeomCam &mgeom1, const Float_t ct1_x, const Float_t ct1_y, const MHillas &hillas2, const MMcEvt &mcevt2, const MGeomCam &mgeom2, const Float_t ct2_x, const Float_t ct2_y)
78{
79 //
80 // Get the direction corresponding to the c.o.g. of the image on
81 // the camera
82 //
83
84 Float_t ct1_cosx_a;
85 Float_t ct1_cosy_a;
86 Float_t ct1_cosz_a; // Direction from ct1 to the shower c.o.g.
87
88 Camera2direction(1e3*mgeom1.GetCameraDist(), mcevt1.GetTelescopePhi(), mcevt1.GetTelescopeTheta(), hillas1.GetMeanX(), hillas1.GetMeanY(), &ct1_cosx_a, &ct1_cosy_a, &ct1_cosz_a);
89
90 //
91 // Now we get another (arbitrary) point along the image long axis,
92 // fMeanX + cosdelta, fMeanY + sindelta, and calculate the direction
93 // to which it corresponds.
94 //
95
96 Float_t ct1_cosx_b;
97 Float_t ct1_cosy_b;
98 Float_t ct1_cosz_b;
99
100 Camera2direction(1e3*mgeom1.GetCameraDist(), mcevt1.GetTelescopePhi(), mcevt1.GetTelescopeTheta(), hillas1.GetMeanX()+hillas1.GetCosDelta(), hillas1.GetMeanY()+hillas1.GetSinDelta(), &ct1_cosx_b, &ct1_cosy_b, &ct1_cosz_b);
101
102 //
103 // The vectorial product of the latter two vectors is a vector
104 // perpendicular to the plane which contains the shower axis and
105 // passes through the telescope center (center of reflector).
106 // The vectorial product of that vector and (0,0,1) is a vector on
107 // the horizontal plane pointing from the telescope center to the
108 // shower core position on the z=0 plane (ground).
109 //
110
111 Float_t ct1_coreVersorX = ct1_cosz_a*ct1_cosx_b - ct1_cosx_a*ct1_cosz_b;
112 Float_t ct1_coreVersorY = ct1_cosz_a*ct1_cosy_b - ct1_cosy_a*ct1_cosz_b;
113
114 //
115 // Now we calculate again the versor, now assuming that the source
116 // direction is paralel to the telescope axis (camera position 0,0)
117 // This increases the precision of the core determination if the showers
118 // actually come from that direction (like for gammas from a point source)
119
120 Camera2direction(1e3*mgeom1.GetCameraDist(), mcevt1.GetTelescopePhi(), mcevt1.GetTelescopeTheta(), 0., 0., &ct1_cosx_b, &ct1_cosy_b, &ct1_cosz_b);
121
122 Float_t ct1_coreVersorX_best = ct1_cosz_a*ct1_cosx_b - ct1_cosx_a*ct1_cosz_b;
123 Float_t ct1_coreVersorY_best = ct1_cosz_a*ct1_cosy_b - ct1_cosy_a*ct1_cosz_b;
124
125 //
126 // Now the second telescope
127 //
128
129 Float_t ct2_cosx_a;
130 Float_t ct2_cosy_a;
131 Float_t ct2_cosz_a; // Direction from ct2 to the shower c.o.g.
132
133
134 Camera2direction(1e3*mgeom2.GetCameraDist(), mcevt2.GetTelescopePhi(), mcevt2.GetTelescopeTheta(), hillas2.GetMeanX(), hillas2.GetMeanY(), &ct2_cosx_a, &ct2_cosy_a, &ct2_cosz_a);
135
136 Float_t ct2_cosx_b;
137 Float_t ct2_cosy_b;
138 Float_t ct2_cosz_b;
139
140 Camera2direction(1e3*mgeom2.GetCameraDist(), mcevt2.GetTelescopePhi(), mcevt2.GetTelescopeTheta(), hillas2.GetMeanX()+hillas2.GetCosDelta(), hillas2.GetMeanY()+hillas2.GetSinDelta(), &ct2_cosx_b, &ct2_cosy_b, &ct2_cosz_b);
141
142
143 Float_t ct2_coreVersorX = ct2_cosz_a*ct2_cosx_b - ct2_cosx_a*ct2_cosz_b;
144 Float_t ct2_coreVersorY = ct2_cosz_a*ct2_cosy_b - ct2_cosy_a*ct2_cosz_b;
145
146
147 Camera2direction(1e3*mgeom2.GetCameraDist(), mcevt2.GetTelescopePhi(), mcevt2.GetTelescopeTheta(), 0., 0., &ct2_cosx_b, &ct2_cosy_b, &ct2_cosz_b);
148
149 Float_t ct2_coreVersorX_best = ct2_cosz_a*ct2_cosx_b - ct2_cosx_a*ct2_cosz_b;
150 Float_t ct2_coreVersorY_best = ct2_cosz_a*ct2_cosy_b - ct2_cosy_a*ct2_cosz_b;
151
152 //
153 // Estimate core position:
154 //
155 Float_t t = ct1_x - ct2_x - ct2_coreVersorX/ct2_coreVersorY*(ct1_y-ct2_y);
156 t /= (ct2_coreVersorX/ct2_coreVersorY*ct1_coreVersorY - ct1_coreVersorX);
157
158 fCoreX = ct1_x + t * ct1_coreVersorX;
159 fCoreY = ct1_y + t * ct1_coreVersorY;
160
161 // fCoreX, fCoreY, fCoreX2, fCoreY2 will have the same units
162 // as ct1_x, ct1_y, ct2_x, ct2_y
163
164
165 //
166 // Now the estimated core position assuming the source is located in
167 // the center of the camera:
168 //
169 t = ct1_x - ct2_x - ct2_coreVersorX_best /
170 ct2_coreVersorY_best*(ct1_y-ct2_y);
171 t /= (ct2_coreVersorX_best/ct2_coreVersorY_best*ct1_coreVersorY_best -
172 ct1_coreVersorX_best);
173
174 fCoreX2 = ct1_x + t * ct1_coreVersorX_best;
175 fCoreY2 = ct1_y + t * ct1_coreVersorY_best;
176
177 //
178 // Be careful, the coordinates in MMcEvt.fCoreX,fCoreY are actually
179 // those of the vector going *from the shower core to the telescope*.
180 // Ours are those of the vector which goes from telescope 1 to the
181 // core estimated core.
182 //
183
184 /////////////////////////////////////////////////////////////////////
185 //
186 // Now estimate the source location on the camera by intersecting
187 // major axis of the ellipses. This assumes both telescopes are
188 // pointing paralel! We introduce the camera scale to account for
189 // the use of telescopes with different focal distances.
190 //
191
192 Float_t scale1 = mgeom1.GetConvMm2Deg();
193 Float_t scale2 = mgeom2.GetConvMm2Deg();
194
195 t = scale2*hillas2.GetMeanY() - scale1*hillas1.GetMeanY() +
196 (scale1*hillas1.GetMeanX() - scale2*hillas2.GetMeanX()) *
197 hillas2.GetSinDelta() / hillas2.GetCosDelta();
198
199 t /= (hillas1.GetSinDelta() -
200 hillas2.GetSinDelta()/hillas2.GetCosDelta()*hillas1.GetCosDelta());
201
202 fSourceX = scale1*hillas1.GetMeanX() + t * hillas1.GetCosDelta();
203 fSourceY = scale1*hillas1.GetMeanY() + t * hillas1.GetSinDelta();
204
205 //
206 // Squared angular distance from reconstructed source position to
207 // camera center.
208 //
209 fTheta2 = fSourceX*fSourceX+fSourceY*fSourceY;
210
211 //
212 // Get the direction corresponding to the intersection of axes
213 //
214
215 Float_t source_direction[3];
216
217 Camera2direction(1e3*mgeom1.GetCameraDist(), mcevt1.GetTelescopePhi(), mcevt1.GetTelescopeTheta(), fSourceX/scale1, fSourceY/scale1, &(source_direction[0]), &(source_direction[1]), &(source_direction[2]));
218
219
220 //
221 // Calculate impact parameters
222 //
223
224 Float_t scalar = (fCoreX-ct1_x)*source_direction[0] +
225 (fCoreY-ct1_y)*source_direction[1];
226 fCT1Impact = sqrt( (fCoreX-ct1_x)*(fCoreX-ct1_x) +
227 (fCoreY-ct1_y)*(fCoreY-ct1_y) -
228 scalar * scalar );
229
230 scalar = (fCoreX-ct2_x)*source_direction[0] +
231 (fCoreY-ct2_y)*source_direction[1];
232 fCT2Impact = sqrt( (fCoreX-ct2_x)*(fCoreX-ct2_x) +
233 (fCoreY-ct2_y)*(fCoreY-ct2_y) -
234 scalar * scalar );
235
236 //
237 // Now calculate i.p. assuming source is point-like and placed in
238 // the center of the camera.
239 //
240 scalar = (fCoreX2-ct1_x)*(-sin(mcevt1.GetTelescopeTheta())*
241 cos(mcevt1.GetTelescopePhi())) +
242 (fCoreY2-ct1_y)*(-sin(mcevt1.GetTelescopeTheta())*
243 sin(mcevt1.GetTelescopePhi()));
244
245 fCT1Impact2 = sqrt( (fCoreX2-ct1_x)*(fCoreX2-ct1_x) +
246 (fCoreY2-ct1_y)*(fCoreY2-ct1_y) -
247 scalar * scalar );
248
249 scalar = (fCoreX2-ct2_x)*(-sin(mcevt2.GetTelescopeTheta())*
250 cos(mcevt2.GetTelescopePhi())) +
251 (fCoreY2-ct2_y)*(-sin(mcevt2.GetTelescopeTheta())*
252 sin(mcevt2.GetTelescopePhi()));
253
254 fCT2Impact2 = sqrt( (fCoreX2-ct2_x)*(fCoreX2-ct2_x) +
255 (fCoreY2-ct2_y)*(fCoreY2-ct2_y) -
256 scalar * scalar );
257
258
259 SetReadyToSave();
260}
261
262// --------------------------------------------------------------------------
263//
264// Transformation of coordinates, from a point on the camera x, y , to
265// the director cosines of the corresponding direction, in the system of
266// coordinates in which X-axis is North, Y-axis is west, and Z-axis
267// points to the zenith. The transformation has been taken from TDAS 01-05,
268// although the final system of coordinates is not the same defined there,
269// but the one defined in Corsika (for convenience).
270//
271// rc is the distance from the reflector center to the camera. CTphi and
272// CTtheta indicate the telescope orientation. The angle CTphi is the
273// azimuth of the vector going along the telescope axis from the camera
274// towards the reflector, measured from the North direction anticlockwise
275// ( being West: phi=pi/2, South: phi=pi, East: phi=3pi/2 )
276//
277// rc and x,y must be given in the same units!
278//
279
280
281void MStereoPar::Camera2direction(Float_t rc, Float_t CTphi, Float_t CTtheta, Float_t x, Float_t y, Float_t* cosx, Float_t* cosy, Float_t* cosz)
282{
283 //
284 // We convert phi to the convention defined in TDAS 01-05
285 //
286 Float_t sinphi = sin(2*TMath::Pi()-CTphi);
287 Float_t cosphi = cos(CTphi);
288 Float_t costheta = cos(CTtheta);
289 Float_t sintheta = sin(CTtheta);
290
291 Float_t xc = x/rc;
292 Float_t yc = y/rc;
293
294 Float_t norm = 1/sqrt(1+xc*xc+yc*yc);
295
296 Float_t xref = xc * norm;
297 Float_t yref = yc * norm;
298 Float_t zref = -1 * norm;
299
300 *cosx = xref * sinphi + yref * costheta*cosphi - zref * sintheta*cosphi;
301 *cosy = -xref * cosphi + yref * costheta*sinphi - zref * sintheta*sinphi;
302 *cosz = yref * sintheta + zref * costheta;
303
304 // Now change from system A of TDAS 01-05 to Corsika system:
305
306 *cosy *= -1;
307 *cosz *= -1;
308
309}
Note: See TracBrowser for help on using the repository browser.