source: trunk/MagicSoft/Mars/mgeom/MGeomCam.cc@ 9313

Last change on this file since 9313 was 9259, checked in by tbretz, 16 years ago
*** empty log message ***
File size: 17.7 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): Thomas Bretz 12/2000 <mailto:tbretz@astro.uni-wuerzburg.de>
19! Author(s): Harald Kornmayer 01/2001
20! Author(s): Markus Gaug 03/2004 <mailto:markus@ifae.es>
21!
22! Copyright: MAGIC Software Development, 2000-2008
23!
24\* ======================================================================== */
25
26/////////////////////////////////////////////////////////////////////////////
27//
28// MGeomCam
29//
30// This is the base class of different camera geometries. It creates
31// a pixel object for a given number of pixels and defines the
32// interface of how to acccess the geometry information.
33//
34// We use a TObjArray for possible future usage (it is much more flexible
35// than a TClonesArray so that it can store more types of pixels (eg
36// fake pixels which are not really existing)
37//
38// Version 1:
39// ----------
40// - first implementation
41//
42// Version 2:
43// ----------
44// - added fPixRatio
45// - added fPixRatioSqrt
46//
47// Version 3:
48// ----------
49// - added fNumAreas
50// - added fNumSectors
51//
52// Version 4:
53// ----------
54// - added fMaxRadius
55// - added fMinRadius
56//
57// Version 5:
58// ----------
59// - added fNumPixInSector
60// - added fNumPixWithAidx
61// - removed fNumSectors
62// - removed fNumAreas
63//
64/////////////////////////////////////////////////////////////////////////////
65#include "MGeomCam.h"
66
67#include <TMath.h> // TMath
68#include <TClass.h> // IsA()->New()
69#include <TArrayI.h> // TArrayI
70#include <TVector2.h> // TVector2
71
72#include "MLog.h"
73#include "MLogManip.h"
74
75#include "MGeomPix.h"
76
77ClassImp(MGeomCam);
78
79using namespace std;
80
81// --------------------------------------------------------------------------
82//
83// Default Constructor. Initializes a Camera Geometry with npix pixels. All
84// pixels are deleted when the corresponding array is deleted.
85//
86MGeomCam::MGeomCam(UInt_t npix, Float_t dist, const char *name, const char *title)
87 : fNumPixels(npix), fCamDist(dist), fConvMm2Deg(kRad2Deg/(dist*1000)),
88 fPixels(npix), fMaxRadius(1), fMinRadius(1), fPixRatio(npix), fPixRatioSqrt(npix)
89{
90 fName = name ? name : "MGeomCam";
91 fTitle = title ? title : "Storage container for a camera geometry";
92
93 //
94 // make sure that the destructor delete all contained objects
95 //
96 fPixels.SetOwner();
97
98 // For root versions <5.18 AddAt is mandatory, for newer
99 // root-version the []-operator can be used safely
100 for (UInt_t i=0; i<npix; i++)
101 fPixels.AddAt(new MGeomPix, i);
102}
103
104// --------------------------------------------------------------------------
105//
106// Copy function. It does intentionally (see workaround in Streamer)
107// not copy name and title.
108//
109void MGeomCam::Copy(TObject &o) const
110{
111 MGeomCam &c = (MGeomCam&)o;
112
113 c.fNumPixels = fNumPixels;
114 c.fCamDist = fCamDist;
115 c.fConvMm2Deg = fConvMm2Deg;
116
117 c.fMaxRadius = fMaxRadius;
118 c.fMinRadius = fMinRadius;
119 c.fPixRatio = fPixRatio;
120 c.fPixRatioSqrt = fPixRatioSqrt;
121
122 c.fNumPixInSector = fNumPixInSector;
123 c.fNumPixWithAidx = fNumPixWithAidx;
124
125 Int_t n = fPixels.GetEntriesFast();
126 Int_t m = c.fPixels.GetEntriesFast();
127
128 c.fPixels.Expand(n);
129
130 for (int i=m; i<n; i++)
131 c.fPixels.AddAt(new MGeomPix, i);
132
133 for (int i=0; i<n; i++)
134 fPixels[i]->Copy(c[i]);
135}
136
137// --------------------------------------------------------------------------
138//
139// Returns a reference of the i-th entry (the pixel with the software idx i)
140// The access is unchecked (for speed reasons!) accesing non existing
141// entries may crash the program!
142//
143MGeomPix &MGeomCam::operator[](Int_t i)
144{
145 return *static_cast<MGeomPix*>(fPixels.UncheckedAt(i));
146}
147
148// --------------------------------------------------------------------------
149//
150// Returns a reference of the i-th entry (the pixel with the software idx i)
151// The access is unchecked (for speed reasons!) accesing non existing
152// entries may crash the program!
153//
154MGeomPix &MGeomCam::operator[](Int_t i) const
155{
156 return *static_cast<MGeomPix*>(fPixels.UncheckedAt(i));
157}
158
159// --------------------------------------------------------------------------
160//
161// Calculate and fill the arrays storing the ratio of the area of a pixel
162// i to the pixel 0 and its square root.
163// The precalculation is done for speed reasons. Having an event the
164// ratio would be calculated at least once for each pixel which is
165// an enormous amount of numerical calculations, which are time
166// consuming and which can be avoided doing the precalculation.
167//
168void MGeomCam::CalcPixRatio()
169{
170 const Double_t a0 = (*this)[0].GetA();
171
172 for (UInt_t i=0; i<fNumPixels; i++)
173 {
174 fPixRatio[i] = a0/(*this)[i].GetA();
175 fPixRatioSqrt[i] = TMath::Sqrt(fPixRatio[i]);
176 }
177}
178
179// --------------------------------------------------------------------------
180//
181// Set the kIsOuterRing flag for all pixels which have a outermost pixel
182// as Next Neighbor and don't have the kIsOutermostRing flag itself.
183//
184void MGeomCam::InitOuterRing()
185{
186 fPixels.R__FOR_EACH(MGeomPix, CheckOuterRing)(*this);
187}
188
189// --------------------------------------------------------------------------
190//
191// Calculate the highest sector index+1 of all pixels, please make sure
192// the the sector numbers are continous.
193//
194void MGeomCam::CalcNumSectors()
195{
196 for (UInt_t i=0; i<fNumPixels; i++)
197 {
198 const Int_t s = (*this)[i].GetSector();
199
200 if (s>=fNumPixInSector.GetSize())
201 fNumPixInSector.Set(s+1);
202
203 fNumPixInSector[s]++;
204 }
205}
206
207// --------------------------------------------------------------------------
208//
209// Calculate the highest area index+1 of all pixels, please make sure
210// the the area indices are continous.
211//
212void MGeomCam::CalcNumAreas()
213{
214 for (UInt_t i=0; i<fNumPixels; i++)
215 {
216 const Int_t s = (*this)[i].GetAidx();
217
218 if (s>=fNumPixWithAidx.GetSize())
219 fNumPixWithAidx.Set(s+1);
220
221 fNumPixWithAidx[s]++;
222 }
223}
224
225// --------------------------------------------------------------------------
226//
227// Calculate the maximum radius of the camera. This is ment for GUI layout.
228//
229void MGeomCam::CalcMaxRadius()
230{
231 fMaxRadius.Set(GetNumAreas()+1);
232 fMinRadius.Set(GetNumAreas()+1);
233
234 for (UInt_t i=0; i<GetNumAreas()+1; i++)
235 {
236 fMaxRadius[i] = 0.;
237 fMinRadius[i] = FLT_MAX;
238 }
239
240 for (UInt_t i=0; i<fNumPixels; i++)
241 {
242 const MGeomPix &pix = (*this)[i];
243
244 const UInt_t s = pix.GetAidx();
245
246 const Float_t d = pix.GetT();
247 const Float_t r = pix.GetDist();
248
249 const Float_t maxr = r + d;
250 const Float_t minr = r>d ? r-d : 0;
251
252 if (maxr>fMaxRadius[s+1])
253 fMaxRadius[s+1] = maxr;
254
255 if (minr<fMinRadius[s+1])
256 fMinRadius[s+1] = minr;
257
258 if (minr<fMinRadius[0])
259 fMinRadius[0] = minr;
260
261 if (maxr>fMaxRadius[0])
262 fMaxRadius[0] = maxr;
263 }
264}
265
266// --------------------------------------------------------------------------
267//
268// sort neighbours from angle of -180 degree to -180 degree
269//
270// where right side of main pixel contains negative degrees
271// and left side positive degrees
272//
273// angle measured from top (0 degree) to bottom (180 degree)
274// ^ - | + //
275// | _ | _ //
276// | / \|/ \ //
277// 30 degree -+- | | | //
278// | / \ / \ / \ //
279// 90 degree -+- | | X | | //
280// | \ / \ / \ / //
281// 150 degree -+- | | | //
282// | \_/|\_/ //
283// | | //
284// | - | + //
285// ------------------> //
286// //
287void MGeomCam::SortNeighbors()
288{
289 for (unsigned int i=0; i<fNumPixels; i++)
290 {
291 MGeomPix &gpix = (*this)[i];
292
293 Double_t phi[6];
294 Int_t idx[7] = { 0, 0, 0, 0, 0, 0, -1 };
295
296 const Int_t n2 = gpix.GetNumNeighbors();
297 for (int j=0; j<n2; j++)
298 {
299 idx[j] = gpix.GetNeighbor(j);
300 phi[j] = (*this)[idx[j]].GetAngle(gpix);
301 }
302
303 Int_t sort[6] = { 6, 6, 6, 6, 6, 6 };
304
305 TMath::Sort(n2, phi, sort, kFALSE);
306
307 gpix.SetNeighbors(idx[sort[0]], idx[sort[1]], idx[sort[2]],
308 idx[sort[3]], idx[sort[4]], idx[sort[5]]);
309 }
310}
311
312// --------------------------------------------------------------------------
313//
314// Returns the distance between the pixels i and j. -1 if an index
315// doesn't exist. The default for j is 0. Assuming that 0 is the index
316// for a central pixel you can get the distance to the camera center.
317//
318Float_t MGeomCam::GetDist(UShort_t i, UShort_t j) const
319{
320 if (i>=fNumPixels || j>=fNumPixels)
321 return -1;
322
323 return (*this)[i].GetDist((*this)[j]);
324}
325
326// --------------------------------------------------------------------------
327//
328// Returns the angle between of pixels i wrt pixel j (default=0). The angle
329// is returned in the range between -pi and pi (atan2) and 2*pi if i or j
330// is out of range.
331//
332Float_t MGeomCam::GetAngle(UShort_t i, UShort_t j) const
333{
334 if (i>=fNumPixels || j>=fNumPixels)
335 return TMath::TwoPi();
336
337 return (*this)[i].GetAngle((*this)[j]);
338}
339
340// --------------------------------------------------------------------------
341//
342// The maximum possible distance from the origin.
343//
344Float_t MGeomCam::GetMaxRadius() const
345{
346 return fMaxRadius[0];
347}
348
349// --------------------------------------------------------------------------
350//
351// The minimum possible distance from the origin.
352//
353Float_t MGeomCam::GetMinRadius() const
354{
355 return fMinRadius[0];
356}
357
358// --------------------------------------------------------------------------
359//
360// Have to call the radii of the subcameras starting to count from 1
361//
362Float_t MGeomCam::GetMaxRadius(const Int_t i) const
363{
364 return i<0 || i>=(Int_t)GetNumAreas() ? -1 : fMaxRadius[i+1];
365}
366
367// --------------------------------------------------------------------------
368//
369// Have to call the radii of the subcameras starting to count from 1
370//
371Float_t MGeomCam::GetMinRadius(const Int_t i) const
372{
373 return i<0 || i>=(Int_t)GetNumAreas() ? -1 : fMinRadius[i+1];
374}
375
376// --------------------------------------------------------------------------
377//
378// returns the ratio of the area of the pixel with index 0 to the pixel
379// with the specified index i. 0 Is returned if the index argument is
380// out of range.
381//
382Float_t MGeomCam::GetPixRatio(UInt_t i) const
383{
384 // Former: (*this)[0].GetA()/(*this)[i].GetA();
385 // The const_cast is necessary to support older root version
386 return i<fNumPixels ? const_cast<TArrayF&>(fPixRatio)[i] : 0;
387}
388
389// --------------------------------------------------------------------------
390//
391// returns the square root of the ratio of the area of the pixel with
392// index 0 to the pixel with the specified index i. 0 Is returned if
393// the index argument is out of range.
394//
395Float_t MGeomCam::GetPixRatioSqrt(UInt_t i) const
396{
397 // The const_cast is necessary to support older root version
398 return i<fNumPixels ? const_cast<TArrayF&>(fPixRatioSqrt)[i] : 0;
399}
400
401// --------------------------------------------------------------------------
402//
403// Check if the position given in the focal plane (so z can be ignored)
404// is a position which might hit the detector. It is meant to be a rough
405// and fast estimate not a precise calculation. All positions dicarded
406// must not hit the detector. All positions accepted might still miss
407// the detector.
408//
409Bool_t MGeomCam::HitDetector(const MQuaternion &v, Double_t offset) const
410{
411 const Double_t max = fMaxRadius[0]/10+offset; // cm --> mm
412 return v.R2()<max*max;
413}
414
415// --------------------------------------------------------------------------
416//
417// Prints the Geometry information of all pixels in the camera.
418// With the option "simple" you can suppress the output of the contents
419// of the individual pixels.
420//
421void MGeomCam::Print(Option_t *o) const
422{
423 //
424 // Print Information about the Geometry of the camera
425 //
426 *fLog << all << " Number of Pixels (" << GetTitle() << "): " << fNumPixels << endl;
427 *fLog << " Number of Sectors: " << GetNumSectors() << " Area-Indices: " << GetNumAreas() << endl;
428 *fLog << " Min.Radius: " << GetMinRadius() << " Max.Radius: " << GetMaxRadius() << endl;
429
430 if (!TString(o).Contains("simple", TString::kIgnoreCase))
431 fPixels.Print();
432}
433
434// --------------------------------------------------------------------------
435//
436// Return the pixel index corresponding to the coordinates given in x, y.
437// The coordinates are given in pixel units (millimeters)
438// If no pixel exists return -1;
439//
440Int_t MGeomCam::GetPixelIdxXY(Float_t x, Float_t y) const
441{
442 for (unsigned int i=0; i<fNumPixels; i++)
443 if ((*this)[i].IsInside(x, y))
444 return i;
445
446 return -1;
447}
448
449Int_t MGeomCam::GetPixelIdx(const TVector2 &v) const
450{
451 return GetPixelIdxXY(v.X(), v.Y());
452}
453
454Int_t MGeomCam::GetPixelIdxDeg(const TVector2 &v) const
455{
456 return GetPixelIdxXYdeg(v.X(), v.Y());
457}
458
459// --------------------------------------------------------------------------
460//
461// Add all indices to arr of all neighbors around pix in a radius r.
462// The center pixel is also returned.
463//
464void MGeomCam::GetNeighbors(TArrayI &arr, const MGeomPix &pix, Float_t r) const
465{
466 arr.Set(GetNumPixels());
467
468 Int_t n = 0;
469
470 for (unsigned int i=0; i<GetNumPixels(); i++)
471 {
472 if (r>TMath::Hypot(pix.GetX()-(*this)[i].GetX(), pix.GetY()-(*this)[i].GetY()))
473 arr[n++] = i;
474 }
475
476 arr.Set(n);
477}
478
479// --------------------------------------------------------------------------
480//
481// Add all indices to arr of all neighbors around idx in a radius r.
482// The center pixel is also returned.
483//
484void MGeomCam::GetNeighbors(TArrayI &arr, UInt_t idx, Float_t r) const
485{
486 if (idx>=GetNumPixels())
487 {
488 arr.Set(0);
489 return;
490 }
491
492 const MGeomPix &pix = (*this)[idx];
493 GetNeighbors(arr, pix, r);
494}
495
496// --------------------------------------------------------------------------
497//
498// Add all pixels to list of all neighbors around pix in a radius r.
499// The center pixel is also returned.
500//
501void MGeomCam::GetNeighbors(TList &arr, const MGeomPix &pix, Float_t r) const
502{
503 for (unsigned int i=0; i<GetNumPixels(); i++)
504 {
505 if (r>TMath::Hypot(pix.GetX()-(*this)[i].GetX(), pix.GetY()-(*this)[i].GetY()))
506 arr.Add(fPixels.UncheckedAt(i));
507 }
508}
509
510// --------------------------------------------------------------------------
511//
512// Add all pixels to list of all neighbors around idx in a radius r.
513// The center pixel is also returned.
514//
515void MGeomCam::GetNeighbors(TList &arr, UInt_t idx, Float_t r) const
516{
517 if (idx>=GetNumPixels())
518 return;
519
520 const MGeomPix &pix = (*this)[idx];
521 GetNeighbors(arr, pix, r);
522}
523
524// --------------------------------------------------------------------------
525//
526// Return direction of p2 w.r.t. p1. For more details
527// see MGeomPix::GetDirection
528//
529Int_t MGeomCam::GetDirection(UInt_t p1, UInt_t p2) const
530{
531 if (p1>fNumPixels || p2>fNumPixels)
532 return -1;
533
534 return operator[](p1).GetDirection(operator[](p2));
535}
536
537// --------------------------------------------------------------------------
538//
539// Get index of neighbor of pixel idx in direction dir, if existing.
540//
541Int_t MGeomCam::GetNeighbor(UInt_t idx, Int_t dir) const
542{
543 if (idx>fNumPixels)
544 return -1;
545
546 const MGeomPix &pix=operator[](idx);
547
548 //
549 // search for the neighbor in the given direction
550 //
551 for (int i=0; i<pix.GetNumNeighbors(); i++)
552 if (GetDirection(idx, pix.GetNeighbor(i))==dir)
553 return pix.GetNeighbor(i);
554
555 return -1;
556}
557
558// --------------------------------------------------------------------------
559//
560// This workaround reproduces (as much as possible) the contents
561// of fPixels which got lost writing MGeomCam in which the
562// fPixels were filles with the []-operator instead of AddAt
563// and a root version previous to 5.18.
564// You try to read broken contents from file if fNumPixels is empty
565// but fNumPixels>0.
566// If you ever read broken contents from a split branch you
567// MUST call this function after reading.
568//
569void MGeomCam::StreamerWorkaround()
570{
571 if (fNumPixels==0 || !fPixels.IsEmpty())
572 return;
573
574 TObject *cam = (TObject*)IsA()->New();
575 cam->Copy(*this);
576 delete cam;
577}
578
579// --------------------------------------------------------------------------
580//
581// This is a custom made streamer. Due to a bug in TObjArray::operator[]
582// old root-versions didn't correctly store the contents of the TObjArray.
583// If such a file is detected (TObjArray empty) a new MGeomCam is created
584// with IsA()->New() and its contents is copied to this. Unfortunately
585// this won't work for all MGeomCam derivatives which need arguments
586// in the constructor and MGeomCam itself (no derivative). Fortunately
587// in prodoction we have never stored anything else than MGeomCamMagic yet.
588// The bug in root can be worked around using AddAt instead of operator[].
589//
590void MGeomCam::Streamer(TBuffer &b)
591{
592 if (b.IsReading())
593 {
594 MGeomCam::Class()->ReadBuffer(b, this);
595 StreamerWorkaround();
596 }
597 else
598 MGeomCam::Class()->WriteBuffer(b, this);
599}
Note: See TracBrowser for help on using the repository browser.