source: trunk/MagicSoft/Mars/mgeom/MGeomPix.cc@ 9385

Last change on this file since 9385 was 9385, checked in by tbretz, 16 years ago
*** empty log message ***
File size: 10.2 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, 12/2000 <mailto:tbretz@astro.uni-wuerzburg.de>
19!
20! Copyright: MAGIC Software Development, 2000-2009
21!
22!
23\* ======================================================================== */
24
25//////////////////////////////////////////////////////////////////////////////
26//
27// MGeomPix
28//
29// This container stores the geometry (position) information of
30// a single pixel together with the information about next neighbors.
31//
32//
33// Version 1:
34// ----------
35// - first implementation
36//
37// Version 2:
38// ----------
39// - added fA
40//
41// Version 3:
42// ----------
43// - added fAidx
44//
45// Version 4:
46// ----------
47// - added fUserBits
48//
49// Version 5:
50// ----------
51// - now derives from MGeom to which some data members have been moved
52//
53////////////////////////////////////////////////////////////////////////////
54#include "MGeomPix.h"
55
56#include <TMath.h>
57#include <TVector2.h>
58#include <TVirtualPad.h>
59
60#include "MLog.h"
61#include "MLogManip.h"
62
63#include "MGeomCam.h"
64#include "MHexagon.h"
65
66ClassImp(MGeomPix);
67
68using namespace std;
69
70const Float_t MGeomPix::gsTan30 = TMath::Sqrt(3)/3; //TMath::Tan(TMath::DegToRad()*30);
71const Float_t MGeomPix::gsTan60 = TMath::Sqrt(3); //TMath::Tan(TMath::DegToRad()*60);
72
73const Float_t MGeomPix::gsCos60 = 0.5; //TMath::Cos(TMath::DegToRad()*60);
74const Float_t MGeomPix::gsSin60 = TMath::Sqrt(3)/2; //TMath::Sin(TMath::DegToRad()*60);
75
76// --------------------------------------------------------------------------
77//
78// Initializes one pixel
79//
80MGeomPix::MGeomPix(Float_t x, Float_t y, Float_t r, UInt_t s, UInt_t a)
81 : MGeom(x, y, s, a), fCosPhi(1), fSinPhi(0)
82{
83 SetD(r);
84 SetNeighbors();
85}
86
87void MGeomPix::SetPhi(Double_t phi)
88{
89 fCosPhi = TMath::Cos(phi);
90 fSinPhi = TMath::Sin(phi);
91}
92
93// --------------------------------------------------------------------------
94//
95// Print the geometry information of one pixel.
96//
97void MGeomPix::Print(Option_t *opt) const
98{
99 MGeom::Print(opt);
100 gLog << " d=" << fD << "mm" << endl;
101}
102
103// ------------------------------------------------------------------------
104//
105// Implementation of PaintPrimitive drwaing a hexagonal pixel
106//
107void MGeomPix::PaintPrimitive(const TAttLine &line, const TAttFill &fill, Double_t scalexy, Double_t scaled) const
108{
109 const_cast<TAttLine&>(line).Modify(); //Change line attributes only if necessary
110 const_cast<TAttFill&>(fill).Modify(); //Change fill area attributes only if necessary
111
112 const Double_t fx = fX*scalexy;
113 const Double_t fy = fY*scalexy;
114 const Double_t fd = fD*scaled;
115
116 const Double_t c = fd*fCosPhi;
117 const Double_t s = fd*fSinPhi;
118
119 const Double_t gsCos60c = gsCos60*c;
120 const Double_t gsCos60s = gsCos60*s;
121
122 const Double_t gsSin60c = gsSin60*c/3;
123 const Double_t gsSin60s = gsSin60*s/3;
124
125 const Double_t gsTan30c = gsTan30*c;
126 const Double_t gsTan30s = gsTan30*s;
127
128
129 Double_t x[7] = {
130 gsCos60c-gsSin60s, -gsTan30s, -gsCos60c-gsSin60s,
131 -gsCos60c+gsSin60s, gsTan30s, gsCos60c+gsSin60s
132 };
133
134 Double_t y[7] = {
135 gsCos60s+gsSin60c, gsTan30c, -gsCos60s+gsSin60c,
136 -gsCos60s-gsSin60c, -gsTan30c, gsCos60s-gsSin60c,
137 };
138
139 for (Int_t i=0; i<6; i++)
140 {
141 // FIXME: Do not rotate fX/fY.
142 // Instead implement MGeomCam::Rotate or similar.
143 x[i] += fx;
144 y[i] += fy;
145 }
146
147 x[6] = x[0];
148 y[6] = y[0];
149
150 //
151 // paint the pixel
152 //
153 if (fill.GetFillColor())
154 gPad->PaintFillArea(6, x, y);
155
156 if (line.GetLineStyle())
157 gPad->PaintPolyLine(7, x, y);
158}
159
160// ------------------------------------------------------------------------
161//
162// compute the distance of a point (px,py) to the Hexagon center in
163// MGeomPix coordinates. Return kTRUE if inside.
164//
165Bool_t MGeomPix::IsInside(Float_t px, Float_t py) const
166{
167
168 //const Double_t disthex = 999999;//TMath::Sqrt(dx*dx + dy*dy);
169
170 const Double_t dh = fD/2;
171
172 //
173 // compute the distance of the Point to the center of the Hexagon
174 //
175 const Double_t dx = px-fX;
176 const Double_t dy = py-fY;
177
178 // FIXME: Derotate
179
180 // DX = dx*cos(-phi) - dy*sin(-phi)
181 // DY = dx*sin(-phi) + dy*cos(-phi)
182
183 const Double_t x = dx*fCosPhi + dy*fSinPhi;
184
185 if (TMath::Abs(x) > dh)
186 return kFALSE; // return disthex
187
188 const Double_t y = -dx*fSinPhi + dy*fCosPhi;
189
190 const Double_t xc = x*gsCos60; // 1/2
191 const Double_t ys = y*gsSin60; // sqrt(3)/2
192
193 //
194 // Now check if point is outside of hexagon; just check x coordinate
195 // in three coordinate systems: the default one, in which two sides of
196 // the hexagon are paralel to the y axis (see camera displays) and two
197 // more, rotated with respect to that one by +- 60 degrees.
198 //
199
200 if (TMath::Abs(xc + ys) > dh)
201 return kFALSE; // return disthex
202
203 if (TMath::Abs(xc - ys) > dh)
204 return kFALSE; // return disthex
205
206 return kTRUE; // return -1
207
208
209}
210
211// ==============================================================================
212/*
213#include <TOrdCollection.h>
214#include "MMath.h"
215#include "../mgui/MHexagon.h" // MHexagon::fgDx
216
217// ------------------------------------------------------------------------
218//
219// Small helper class to allow fast sorting of TVector2 by angle
220//
221class HVector2 : public TVector2
222{
223 Double_t fAngle;
224public:
225 HVector2() : TVector2() { }
226 HVector2(const TVector2 &v) : TVector2(v) { }
227 HVector2(Double_t x, Double_t y) : TVector2 (x, y) { }
228 void CalcAngle() { fAngle = Phi(); }
229 Bool_t IsSortable() const { return kTRUE; }
230 Int_t Compare(const TObject *obj) const { return ((HVector2*)obj)->fAngle>fAngle ? 1 : -1; }
231 Double_t GetAngle() const { return fAngle; }
232};
233
234// ------------------------------------------------------------------------
235//
236// Calculate the edge points of the intersection area of two hexagons.
237// The points are added as TVector2 to the TOrdCollection.
238// The user is responsible of delete the objects.
239//
240void MGeomPix::GetIntersectionBorder(TOrdCollection &col, const MGeomPix &hex) const
241{
242// if (fPhi!=0)
243// {
244// gLog << warn << "MGeomPix::GetIntersectionBorder: Only phi=0 supported yet." << endl;
245// return;
246// }
247
248 Bool_t inside0[6], inside1[6];
249
250 HVector2 v0[6];
251 HVector2 v1[6];
252
253 Int_t cnt0=0;
254 Int_t cnt1=0;
255
256 // Calculate teh edges of each hexagon and whether this edge
257 // is inside the other hexgon or not
258 for (int i=0; i<6; i++)
259 {
260 const Double_t x0 = fX+MHexagon::fgDx[i]*fD;
261 const Double_t y0 = fY+MHexagon::fgDy[i]*fD;
262
263 const Double_t x1 = hex.fX+MHexagon::fgDx[i]*hex.fD;
264 const Double_t y1 = hex.fY+MHexagon::fgDy[i]*hex.fD;
265
266 v0[i].Set(x0, y0);
267 v1[i].Set(x1, y1);
268
269 inside0[i] = hex.DistanceToPrimitive(x0, y0) < 0;
270 inside1[i] = DistanceToPrimitive(x1, y1) < 0;
271
272 if (inside0[i])
273 {
274 col.Add(new HVector2(v0[i]));
275 cnt0++;
276 }
277 if (inside1[i])
278 {
279 col.Add(new HVector2(v1[i]));
280 cnt1++;
281 }
282 }
283
284 if (cnt0==0 || cnt1==0)
285 return;
286
287 // No calculate which vorder lines intersect
288 Bool_t iscross0[6], iscross1[6];
289 for (int i=0; i<6; i++)
290 {
291 iscross0[i] = (inside0[i] && !inside0[(i+1)%6]) || (!inside0[i] && inside0[(i+1)%6]);
292 iscross1[i] = (inside1[i] && !inside1[(i+1)%6]) || (!inside1[i] && inside1[(i+1)%6]);
293 }
294
295 // Calculate the border points of our intersection area
296 for (int i=0; i<6; i++)
297 {
298 // Skip non intersecting lines
299 if (!iscross0[i])
300 continue;
301
302 for (int j=0; j<6; j++)
303 {
304 // Skip non intersecting lines
305 if (!iscross1[j])
306 continue;
307
308 const TVector2 p = MMath::GetIntersectionPoint(v0[i], v0[(i+1)%6], v1[j], v1[(j+1)%6]);
309 if (hex.DistanceToPrimitive(p.X(), p.Y())<1e-9)
310 col.Add(new HVector2(p));
311 }
312 }
313}
314
315// ------------------------------------------------------------------------
316//
317// Calculate the overlapping area of the two hexagons.
318//
319Double_t MGeomPix::CalcOverlapArea(const MGeomPix &cam) const
320{
321// if (fPhi!=0)
322// {
323// gLog << warn << "MGeomPix::CalcOverlapArea: Only phi=0 supported yet." << endl;
324// return -1;
325// }
326
327 TOrdCollection col;
328 col.SetOwner();
329
330 GetIntersectionBorder(col, cam);
331
332 // Check if there is an intersection to proceed with
333 const Int_t n = col.GetEntries();
334 if (n==0)
335 return 0;
336
337 // Calculate the center of gravity
338 TVector2 cog;
339
340 TIter Next(&col);
341 HVector2 *v=0;
342 while ((v=(HVector2*)Next()))
343 cog += *v;
344 cog /= n;
345
346 // Shift the figure to its center-og-gravity and
347 // calculate the angle of the connection line between the
348 // border points of our intersesction area and its cog
349 Next.Reset();
350 while ((v=(HVector2*)Next()))
351 {
352 *v -= cog;
353 v->CalcAngle();
354 }
355
356 // Sort these points by this angle
357 col.Sort();
358
359 // Now sum up the area of all the triangles between two
360 // following points and the cog.
361 Double_t A = 0;
362 for (int i=0; i<n; i++)
363 {
364 // Vectors from cog to two nearby border points
365 const HVector2 *v1 = (HVector2*)col.At(i);
366 const HVector2 *v2 = (HVector2*)col.At((i+1)%n);
367
368 // Angle between both vectors
369 const Double_t a = fmod(v1->GetAngle()-v2->GetAngle()+TMath::TwoPi(), TMath::TwoPi());
370
371 // Length of both vectors
372 const Double_t d1 = v1->Mod();
373 const Double_t d2 = v2->Mod();
374
375 A += d1*d2/2*sin(a);
376 }
377 return A;
378}
379*/
Note: See TracBrowser for help on using the repository browser.