source: trunk/Mars/mgeom/MGeomCam.cc@ 10090

Last change on this file since 10090 was 10090, checked in by tbretz, 12 years ago
Implemented improved handling for the obscuration by the detector.
File size: 20.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 "MGeom.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)), fPixels(npix),
88 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
99
100// --------------------------------------------------------------------------
101//
102// To make sure that everything is correctly initialized when
103// cloned we call InitGeometry after cloning. This repairs
104// a few I/O problems mainly with files written with older
105// MARS versions (see also StreamerWorkaround)
106//
107TObject *MGeomCam::Clone(const char *newname) const
108{
109 MGeomCam *clone = (MGeomCam*)MParContainer::Clone(newname);
110 clone->InitGeometry();
111 return clone;
112}
113
114// --------------------------------------------------------------------------
115//
116// Copy function. It does intentionally (see workaround in Streamer)
117// not copy name and title.
118//
119void MGeomCam::Copy(TObject &o) const
120{
121 MGeomCam &c = (MGeomCam&)o;
122
123 c.fNumPixels = fNumPixels;
124 c.fCamDist = fCamDist;
125 c.fConvMm2Deg = fConvMm2Deg;
126
127 c.fMaxRadius = fMaxRadius;
128 c.fMinRadius = fMinRadius;
129 c.fPixRatio = fPixRatio;
130 c.fPixRatioSqrt = fPixRatioSqrt;
131
132 c.fNumPixInSector = fNumPixInSector;
133 c.fNumPixWithAidx = fNumPixWithAidx;
134
135 Int_t n = fPixels.GetEntriesFast();
136 Int_t m = c.fPixels.GetEntriesFast();
137
138 c.fPixels.Delete();
139 c.fPixels.Expand(n);
140
141 for (int i=m; i<n; i++)
142 c.fPixels.AddAt(fPixels[i]->Clone(), i);
143}
144
145// --------------------------------------------------------------------------
146//
147// Returns a reference of the i-th entry (the pixel with the software idx i)
148// The access is unchecked (for speed reasons!) accesing non existing
149// entries may crash the program!
150//
151MGeom &MGeomCam::operator[](Int_t i)
152{
153 return *static_cast<MGeom*>(fPixels.UncheckedAt(i));
154}
155
156// --------------------------------------------------------------------------
157//
158// Returns a reference of the i-th entry (the pixel with the software idx i)
159// The access is unchecked (for speed reasons!) accesing non existing
160// entries may crash the program!
161//
162MGeom &MGeomCam::operator[](Int_t i) const
163{
164 return *static_cast<MGeom*>(fPixels.UncheckedAt(i));
165}
166
167// --------------------------------------------------------------------------
168//
169// Calculate and fill the arrays storing the ratio of the area of a pixel
170// i to the pixel 0 and its square root.
171// The precalculation is done for speed reasons. Having an event the
172// ratio would be calculated at least once for each pixel which is
173// an enormous amount of numerical calculations, which are time
174// consuming and which can be avoided doing the precalculation.
175//
176void MGeomCam::CalcPixRatio()
177{
178 const Double_t a0 = (*this)[0].GetA();
179
180 for (UInt_t i=0; i<fNumPixels; i++)
181 {
182 fPixRatio[i] = a0/(*this)[i].GetA();
183 fPixRatioSqrt[i] = TMath::Sqrt(fPixRatio[i]);
184 }
185}
186
187// --------------------------------------------------------------------------
188//
189// Set the kIsOuterRing flag for all pixels which have a outermost pixel
190// as Next Neighbor and don't have the kIsOutermostRing flag itself.
191//
192void MGeomCam::InitOuterRing()
193{
194 fPixels.R__FOR_EACH(MGeom, CheckOuterRing)(*this);
195}
196
197// --------------------------------------------------------------------------
198//
199// Calculate the highest sector index+1 of all pixels, please make sure
200// the the sector numbers are continous.
201//
202void MGeomCam::CalcNumSectors()
203{
204 for (UInt_t i=0; i<fNumPixels; i++)
205 {
206 const Int_t s = (*this)[i].GetSector();
207
208 if (s>=fNumPixInSector.GetSize())
209 fNumPixInSector.Set(s+1);
210
211 fNumPixInSector[s]++;
212 }
213}
214
215// --------------------------------------------------------------------------
216//
217// Calculate the highest area index+1 of all pixels, please make sure
218// the the area indices are continous.
219//
220void MGeomCam::CalcNumAreas()
221{
222 for (UInt_t i=0; i<fNumPixels; i++)
223 {
224 const Int_t s = (*this)[i].GetAidx();
225
226 if (s>=fNumPixWithAidx.GetSize())
227 fNumPixWithAidx.Set(s+1);
228
229 fNumPixWithAidx[s]++;
230 }
231}
232
233// --------------------------------------------------------------------------
234//
235// Calculate the maximum radius of the camera. This is ment for GUI layout.
236//
237void MGeomCam::CalcMaxRadius()
238{
239 fMaxRadius.Set(GetNumAreas()+1);
240 fMinRadius.Set(GetNumAreas()+1);
241
242 for (UInt_t i=0; i<GetNumAreas()+1; i++)
243 {
244 fMaxRadius[i] = 0.;
245 fMinRadius[i] = FLT_MAX;
246 }
247
248 for (UInt_t i=0; i<fNumPixels; i++)
249 {
250 const MGeom &pix = (*this)[i];
251
252 const UInt_t s = pix.GetAidx();
253
254 const Float_t d = pix.GetT();
255 const Float_t r = pix.GetDist();
256
257 const Float_t maxr = r + d;
258 const Float_t minr = r>d ? r-d : 0;
259
260 if (maxr>fMaxRadius[s+1])
261 fMaxRadius[s+1] = maxr;
262
263 if (minr<fMinRadius[s+1])
264 fMinRadius[s+1] = minr;
265
266 if (minr<fMinRadius[0])
267 fMinRadius[0] = minr;
268
269 if (maxr>fMaxRadius[0])
270 fMaxRadius[0] = maxr;
271 }
272}
273
274// --------------------------------------------------------------------------
275//
276// sort neighbours from angle of -180 degree to -180 degree
277//
278// where right side of main pixel contains negative degrees
279// and left side positive degrees
280//
281// angle measured from top (0 degree) to bottom (180 degree)
282// ^ - | + //
283// | _ | _ //
284// | / \|/ \ //
285// 30 degree -+- | | | //
286// | / \ / \ / \ //
287// 90 degree -+- | | X | | //
288// | \ / \ / \ / //
289// 150 degree -+- | | | //
290// | \_/|\_/ //
291// | | //
292// | - | + //
293// ------------------> //
294// //
295void MGeomCam::SortNeighbors()
296{
297 for (unsigned int i=0; i<fNumPixels; i++)
298 {
299 MGeom &gpix = (*this)[i];
300
301 Double_t phi[6];
302 Int_t idx[7] = { 0, 0, 0, 0, 0, 0, -1 };
303
304 const Int_t n2 = gpix.GetNumNeighbors();
305 for (int j=0; j<n2; j++)
306 {
307 idx[j] = gpix.GetNeighbor(j);
308 phi[j] = (*this)[idx[j]].GetAngle(gpix);
309 }
310
311 Int_t sort[6] = { 6, 6, 6, 6, 6, 6 };
312
313 TMath::Sort(n2, phi, sort, kFALSE);
314
315 gpix.SetNeighbors(idx[sort[0]], idx[sort[1]], idx[sort[2]],
316 idx[sort[3]], idx[sort[4]], idx[sort[5]]);
317 }
318}
319
320// --------------------------------------------------------------------------
321//
322// Return the total area of the camera
323//
324Float_t MGeomCam::GetA() const
325{
326 Double_t A = 0;
327 for (unsigned int i=0; i<fNumPixels; i++)
328 A += (*this)[i].GetA();
329
330 return A;
331}
332
333// --------------------------------------------------------------------------
334//
335// Returns the distance between the pixels i and j. -1 if an index
336// doesn't exist. The default for j is 0. Assuming that 0 is the index
337// for a central pixel you can get the distance to the camera center.
338//
339Float_t MGeomCam::GetDist(UShort_t i, UShort_t j) const
340{
341 if (i>=fNumPixels || j>=fNumPixels)
342 return -1;
343
344 return (*this)[i].GetDist((*this)[j]);
345}
346
347// --------------------------------------------------------------------------
348//
349// Returns the angle between of pixels i wrt pixel j (default=0). The angle
350// is returned in the range between -pi and pi (atan2) and 2*pi if i or j
351// is out of range.
352//
353Float_t MGeomCam::GetAngle(UShort_t i, UShort_t j) const
354{
355 if (i>=fNumPixels || j>=fNumPixels)
356 return TMath::TwoPi();
357
358 return (*this)[i].GetAngle((*this)[j]);
359}
360
361// --------------------------------------------------------------------------
362//
363// The maximum possible distance from the origin.
364//
365Float_t MGeomCam::GetMaxRadius() const
366{
367 return fMaxRadius[0];
368}
369
370// --------------------------------------------------------------------------
371//
372// The minimum possible distance from the origin.
373//
374Float_t MGeomCam::GetMinRadius() const
375{
376 return fMinRadius[0];
377}
378
379// --------------------------------------------------------------------------
380//
381// Have to call the radii of the subcameras starting to count from 1
382//
383Float_t MGeomCam::GetMaxRadius(const Int_t i) const
384{
385 return i<0 || i>=(Int_t)GetNumAreas() ? -1 : fMaxRadius[i+1];
386}
387
388// --------------------------------------------------------------------------
389//
390// Have to call the radii of the subcameras starting to count from 1
391//
392Float_t MGeomCam::GetMinRadius(const Int_t i) const
393{
394 return i<0 || i>=(Int_t)GetNumAreas() ? -1 : fMinRadius[i+1];
395}
396
397// --------------------------------------------------------------------------
398//
399// returns the ratio of the area of the pixel with index 0 to the pixel
400// with the specified index i. 0 Is returned if the index argument is
401// out of range.
402//
403Float_t MGeomCam::GetPixRatio(UInt_t i) const
404{
405 // Former: (*this)[0].GetA()/(*this)[i].GetA();
406 // The const_cast is necessary to support older root version
407 return i<fNumPixels ? const_cast<TArrayF&>(fPixRatio)[i] : 0;
408}
409
410// --------------------------------------------------------------------------
411//
412// returns the square root of the ratio of the area of the pixel with
413// index 0 to the pixel with the specified index i. 0 Is returned if
414// the index argument is out of range.
415//
416Float_t MGeomCam::GetPixRatioSqrt(UInt_t i) const
417{
418 // The const_cast is necessary to support older root version
419 return i<fNumPixels ? const_cast<TArrayF&>(fPixRatioSqrt)[i] : 0;
420}
421
422// --------------------------------------------------------------------------
423//
424// Check if the photon which is flying along the trajectory u has passed
425// (or will pass) the frame of the camera (and consequently get
426// absorbed). The position p and direction u must be in the
427// telescope coordinate frame, which is z parallel to the focal plane,
428// x to the right and y upwards, looking from the mirror towards the camera.
429//
430// The units are cm.
431//
432Bool_t MGeomCam::HitFrame(MQuaternion p, const MQuaternion &u, Double_t margin) const
433{
434 if (margin<=0)
435 return HitFrame(p, u);
436
437 // z is defined from the mirror (0) to the camera (z>0).
438 // Thus we just propagate to the focal plane (z=fDist)
439 p.PropagateZ(u, GetCameraDist()*100); // m->cm
440 return p.R()<margin;
441}
442
443// --------------------------------------------------------------------------
444//
445// Check if the position given in the focal plane (so z can be ignored)
446// is a position which might hit the detector. It is meant to be a rough
447// and fast estimate not a precise calculation. All positions dicarded
448// must not hit the detector. All positions accepted might still miss
449// the detector.
450//
451Bool_t MGeomCam::HitDetector(const MQuaternion &v, Double_t offset) const
452{
453 const Double_t max = fMaxRadius[0]/10+offset; // cm --> mm
454 return v.R2()<max*max;
455}
456
457// --------------------------------------------------------------------------
458//
459// Prints the Geometry information of all pixels in the camera.
460// With the option "simple" you can suppress the output of the contents
461// of the individual pixels.
462//
463void MGeomCam::Print(Option_t *o) const
464{
465 //
466 // Print Information about the Geometry of the camera
467 //
468 *fLog << all << " Number of Pixels (" << GetTitle() << "): " << fNumPixels << endl;
469 *fLog << " Number of Sectors: " << GetNumSectors() << " Area-Indices: " << GetNumAreas() << endl;
470 *fLog << " Min.Radius: " << GetMinRadius() << " Max.Radius: " << GetMaxRadius() << endl;
471
472 if (!TString(o).Contains("simple", TString::kIgnoreCase))
473 fPixels.Print();
474}
475
476// --------------------------------------------------------------------------
477//
478// Return the pixel index corresponding to the coordinates given in x, y.
479// The coordinates are given in pixel units (millimeters)
480// If no pixel exists return -1;
481//
482Int_t MGeomCam::GetPixelIdxXY(Float_t x, Float_t y) const
483{
484 for (unsigned int i=0; i<fNumPixels; i++)
485 if ((*this)[i].IsInside(x, y))
486 return i;
487
488 return -1;
489}
490
491Int_t MGeomCam::GetPixelIdx(const TVector2 &v) const
492{
493 return GetPixelIdxXY(v.X(), v.Y());
494}
495
496Int_t MGeomCam::GetPixelIdxDeg(const TVector2 &v) const
497{
498 return GetPixelIdxXYdeg(v.X(), v.Y());
499}
500
501// --------------------------------------------------------------------------
502//
503// Add all indices to arr of all neighbors around pix in a radius r.
504// The center pixel is also returned.
505//
506void MGeomCam::GetNeighbors(TArrayI &arr, const MGeom &pix, Float_t r) const
507{
508 arr.Set(GetNumPixels());
509
510 Int_t n = 0;
511
512 for (unsigned int i=0; i<GetNumPixels(); i++)
513 {
514 if (r>TMath::Hypot(pix.GetX()-(*this)[i].GetX(), pix.GetY()-(*this)[i].GetY()))
515 arr[n++] = i;
516 }
517
518 arr.Set(n);
519}
520
521// --------------------------------------------------------------------------
522//
523// Add all indices to arr of all neighbors around idx in a radius r.
524// The center pixel is also returned.
525//
526void MGeomCam::GetNeighbors(TArrayI &arr, UInt_t idx, Float_t r) const
527{
528 if (idx>=GetNumPixels())
529 {
530 arr.Set(0);
531 return;
532 }
533
534 const MGeom &pix = (*this)[idx];
535 GetNeighbors(arr, pix, r);
536}
537
538// --------------------------------------------------------------------------
539//
540// Add all pixels to list of all neighbors around pix in a radius r.
541// The center pixel is also returned.
542//
543void MGeomCam::GetNeighbors(TList &arr, const MGeom &pix, Float_t r) const
544{
545 for (unsigned int i=0; i<GetNumPixels(); i++)
546 {
547 if (r>TMath::Hypot(pix.GetX()-(*this)[i].GetX(), pix.GetY()-(*this)[i].GetY()))
548 arr.Add(fPixels.UncheckedAt(i));
549 }
550}
551
552// --------------------------------------------------------------------------
553//
554// Add all pixels to list of all neighbors around idx in a radius r.
555// The center pixel is also returned.
556//
557void MGeomCam::GetNeighbors(TList &arr, UInt_t idx, Float_t r) const
558{
559 if (idx>=GetNumPixels())
560 return;
561
562 const MGeom &pix = (*this)[idx];
563 GetNeighbors(arr, pix, r);
564}
565
566// --------------------------------------------------------------------------
567//
568// Return direction of p2 w.r.t. p1. For more details
569// see MGeom::GetDirection
570//
571Int_t MGeomCam::GetDirection(UInt_t p1, UInt_t p2) const
572{
573 if (p1>fNumPixels || p2>fNumPixels)
574 return -1;
575
576 return operator[](p1).GetDirection(operator[](p2));
577}
578
579// --------------------------------------------------------------------------
580//
581// Get index of neighbor of pixel idx in direction dir, if existing.
582//
583Int_t MGeomCam::GetNeighbor(UInt_t idx, Int_t dir) const
584{
585 if (idx>fNumPixels)
586 return -1;
587
588 const MGeom &pix=operator[](idx);
589
590 //
591 // search for the neighbor in the given direction
592 //
593 for (int i=0; i<pix.GetNumNeighbors(); i++)
594 if (GetDirection(idx, pix.GetNeighbor(i))==dir)
595 return pix.GetNeighbor(i);
596
597 return -1;
598}
599
600// --------------------------------------------------------------------------
601//
602// This fills the next neighbor information from a table into the pixel
603// objects. An overloaded function must call SortNeighbors()
604// at the end.
605//
606// All pixels with a center withing 1.75*GetT() are considered neighbors,
607// Only six neighbors are allowed.
608//
609void MGeomCam::CreateNN()
610{
611 TArrayI nn(6);
612
613 for (UInt_t i=0; i<GetNumPixels(); i++)
614 {
615 MGeom &pix = (*this)[i];
616
617 Int_t k = 0;
618 nn.Reset(-1);
619
620 for (UInt_t j=0; j<GetNumPixels(); j++)
621 {
622 if (i==j)
623 continue;
624
625 if (pix.GetDist((*this)[j])>pix.GetT()*1.75)
626 continue;
627
628 if (k==6)
629 {
630 *fLog << err << "ERROR - MGeomCam::CreateNN: Pixel " << j << " has too many neighbors." << endl;
631 break;
632 }
633
634 nn[k++] = j;
635 }
636
637 pix.SetNeighbors(nn[0], nn[1], nn[2], nn[3], nn[4], nn[5]);
638 }
639
640 SortNeighbors();
641}
642
643// --------------------------------------------------------------------------
644//
645// This workaround reproduces (as much as possible) the contents
646// of fPixels which got lost writing MGeomCam in which the
647// fPixels were filles with the []-operator instead of AddAt
648// and a root version previous to 5.18.
649// You try to read broken contents from file if fNumPixels is empty
650// but fNumPixels>0.
651// If you ever read broken contents from a split branch you
652// MUST call this function after reading.
653//
654// Furthermore since we moved all contents of MGeomPix to its
655// new base class the contents of fNeighbors is not read anymore
656// although all other contents is read correctly. Therefore we have
657// top recreate the neighbor table.
658//
659void MGeomCam::StreamerWorkaround()
660{
661 if (fNumPixels==0)
662 return;
663
664
665 if (!fPixels.IsEmpty())
666 {
667 CreateNN();
668 return;
669 }
670
671 const TObject *cam = (TObject*)IsA()->New();
672 cam->Copy(*this);
673 delete cam;
674}
675
676// --------------------------------------------------------------------------
677//
678// This is a custom made streamer. Due to a bug in TObjArray::operator[]
679// old root-versions didn't correctly store the contents of the TObjArray.
680// If such a file is detected (TObjArray empty) a new MGeomCam is created
681// with IsA()->New() and its contents is copied to this. Unfortunately
682// this won't work for all MGeomCam derivatives which need arguments
683// in the constructor and MGeomCam itself (no derivative). Fortunately
684// in prodoction we have never stored anything else than MGeomCamMagic yet.
685// The bug in root can be worked around using AddAt instead of operator[].
686//
687void MGeomCam::Streamer(TBuffer &b)
688{
689 if (b.IsReading())
690 {
691 MGeomCam::Class()->ReadBuffer(b, this);
692 StreamerWorkaround();
693 }
694 else
695 MGeomCam::Class()->WriteBuffer(b, this);
696}
697
698// --------------------------------------------------------------------------
699//
700// This deletes the pixel i and sets it to a clone of pix.
701// If i>=fNumPixels it is just ignored.
702//
703void MGeomCam::SetAt(UInt_t i, const MGeom &pix)
704{
705 if (i>=fNumPixels)
706 return;
707
708 if (fPixels[i])
709 delete fPixels.RemoveAt(i);
710
711 // For root versions <5.18 AddAt is mandatory, for newer
712 // root-version the []-operator can be used safely
713 fPixels.AddAt(pix.Clone(), i);
714}
Note: See TracBrowser for help on using the repository browser.