source: branches/AddingGoogleTestEnvironment/msignal/MSignalCam.cc@ 19690

Last change on this file since 19690 was 9573, checked in by tbretz, 15 years ago
*** empty log message ***
File size: 18.3 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, 03/2005 <mailto:tbretz@astro.uni-wuerzburg.de>
19!
20! Copyright: MAGIC Software Development, 2000-2010
21!
22! Version 2:
23! ----------
24!
25! Made persistent (the correct contets is not guranteed, because it is
26! calculated from the outside (e.g. in MImgCleanStd))
27!
28! Short_t fNumIslands; //!
29! Short_t fNumSinglePixels; //!
30! Float_t fSizeSinglePixels; //!
31! Float_t fSizeSubIslands; //!
32! Float_t fSizeMainIsland; //!
33!
34!
35\* ======================================================================== */
36
37/////////////////////////////////////////////////////////////////////////////
38//
39// MSignalCam
40//
41// Class Version 1:
42// ----------------
43// - first version
44//
45/////////////////////////////////////////////////////////////////////////////
46#include "MSignalCam.h"
47
48#include <math.h>
49#include <limits.h>
50#include <fstream>
51
52#include <TArrayI.h>
53#include <TArrayD.h>
54
55#include "MLog.h"
56#include "MLogManip.h"
57
58#include "MGeomCam.h"
59#include "MGeom.h"
60
61ClassImp(MSignalCam);
62ClassImp(MSignalCamIter);
63
64using namespace std;
65
66// --------------------------------------------------------------------------
67//
68// Creates a MSignalPix object for each pixel in the event
69//
70MSignalCam::MSignalCam(const char *name, const char *title) /*: fNumPixels(0)*/
71{
72 fName = name ? name : "MSignalCam";
73 fTitle = title ? title : "(Number of Photon)-Event Information";
74
75 fPixels = new TClonesArray("MSignalPix", 0);
76
77 Reset();
78}
79
80// --------------------------------------------------------------------------
81//
82// Copy contants of obj into this instance.
83//
84void MSignalCam::Copy(TObject &obj) const
85{
86 MSignalCam &cam = static_cast<MSignalCam&>(obj);
87
88 cam.fNumIslands = fNumIslands;
89 cam.fNumSinglePixels = fNumSinglePixels;
90 cam.fSizeSinglePixels = fSizeSinglePixels;
91 cam.fSizeSubIslands = fSizeSubIslands;
92 cam.fSizeMainIsland = fSizeMainIsland;
93
94 cam.fNumPixelsSaturatedHiGain = fNumPixelsSaturatedHiGain;
95 cam.fNumPixelsSaturatedLoGain = fNumPixelsSaturatedLoGain;
96
97 cam.fPixels->Delete();
98
99 const UInt_t n = GetNumPixels();
100
101 cam.InitSize(n);
102
103 for (UInt_t i=0; i<n; i++)
104 new ((*cam.fPixels)[i]) MSignalPix((*this)[i]);
105}
106
107// --------------------------------------------------------------------------
108//
109// reset counter and delete netries in list.
110//
111void MSignalCam::Reset()
112{
113 fNumSinglePixels = 0;
114 fSizeSinglePixels = 0;
115 fSizeSubIslands = 0;
116 fSizeMainIsland = 0;
117 fNumIslands = -1;
118
119 fNumPixelsSaturatedHiGain = -1;
120 fNumPixelsSaturatedLoGain = -1;
121
122 fPixels->R__FOR_EACH(TObject, Clear)();
123}
124
125// --------------------------------------------------------------------------
126//
127// Dump the cerenkov photon event to *fLog
128//
129void MSignalCam::Print(Option_t *) const
130{
131 const Int_t entries = fPixels->GetEntries();
132
133 *fLog << GetDescriptor() << dec << endl;
134 *fLog << " Number of Pixels: " << GetNumPixels() << "(" << entries << ")" << endl;
135
136 for (Int_t i=0; i<entries; i++ )
137 (*this)[i].Print();
138}
139
140// --------------------------------------------------------------------------
141//
142// Count and return the number of unmapped pixels
143//
144Int_t MSignalCam::GetNumPixelsUnmapped() const
145{
146 const UInt_t n = GetNumPixels();
147
148 if (n <= 0)
149 return -1;
150
151 Int_t cnt=0;
152 for (UInt_t i=0; i<n; i++)
153 {
154 if ((*this)[i].IsPixelUnmapped())
155 cnt++;
156 }
157
158 return cnt;
159}
160
161// --------------------------------------------------------------------------
162//
163// get the minimum number of photons of all valid pixels in the list
164// If you specify a geometry the number of photons is weighted with the
165// area of the pixel
166//
167Float_t MSignalCam::GetNumPhotonsMin(const MGeomCam *geom) const
168{
169 const UInt_t n = GetNumPixels();
170
171 if (n <= 0)
172 return -5.;
173
174 Float_t minval = FLT_MAX;
175
176 for (UInt_t i=0; i<n; i++)
177 {
178 const MSignalPix &pix = (*this)[i];
179 if (!pix.IsPixelUsed())
180 continue;
181
182 Float_t testval = pix.GetNumPhotons();
183
184 if (geom)
185 testval *= geom->GetPixRatio(i);
186
187 if (testval < minval)
188 minval = testval;
189 }
190
191 return minval;
192}
193
194// --------------------------------------------------------------------------
195//
196// get the maximum number of photons of all valid pixels in the list
197// If you specify a geometry the number of photons is weighted with the
198// area of the pixel
199//
200Float_t MSignalCam::GetNumPhotonsMax(const MGeomCam *geom) const
201{
202 const UInt_t n = GetNumPixels();
203
204 if (n <= 0)
205 return 50.;
206
207 Float_t maxval = -FLT_MAX;
208
209 for (UInt_t i=0; i<n; i++)
210 {
211 const MSignalPix &pix = (*this)[i];
212 if (!pix.IsPixelUsed())
213 continue;
214
215 Float_t testval = pix.GetNumPhotons();
216 if (geom)
217 testval *= geom->GetPixRatio(i);
218
219 if (testval > maxval)
220 maxval = testval;
221 }
222 return maxval;
223}
224
225// --------------------------------------------------------------------------
226//
227// get the minimum ratio of photons/error
228//
229Float_t MSignalCam::GetRatioMin(const MGeomCam *geom) const
230{
231 const UInt_t n = GetNumPixels();
232
233 if (n <= 0)
234 return -5.;
235
236 Float_t minval = FLT_MAX;
237
238 for (UInt_t i=0; i<n; i++)
239 {
240 const MSignalPix &pix = (*this)[i];
241 if (!pix.IsPixelUsed())
242 continue;
243
244 Float_t testval = pix.GetNumPhotons()/pix.GetErrorPhot();
245 if (geom)
246 testval *= geom->GetPixRatioSqrt(i);
247
248 if (testval < minval)
249 minval = testval;
250 }
251
252 return minval;
253}
254
255// --------------------------------------------------------------------------
256//
257// get the maximum ratio of photons/error
258//
259Float_t MSignalCam::GetRatioMax(const MGeomCam *geom) const
260{
261 const UInt_t n = GetNumPixels();
262
263 if (n <= 0)
264 return -5.;
265
266 Float_t maxval = -FLT_MAX;
267
268 for (UInt_t i=0; i<n; i++)
269 {
270 const MSignalPix &pix = (*this)[i];
271 if (!pix.IsPixelUsed())
272 continue;
273
274 Float_t testval = pix.GetNumPhotons()/pix.GetErrorPhot();
275 if (geom)
276 testval *= geom->GetPixRatioSqrt(i);
277
278 if (testval > maxval)
279 maxval = testval;
280 }
281
282 return maxval;
283}
284
285// --------------------------------------------------------------------------
286//
287// get the minimum of error
288// If you specify a geometry the number of photons is weighted with the
289// area of the pixel
290//
291Float_t MSignalCam::GetErrorPhotMin(const MGeomCam *geom) const
292{
293 const UInt_t n = GetNumPixels();
294
295 if (n <= 0)
296 return 50.;
297
298 Float_t minval = FLT_MAX;
299
300 for (UInt_t i=0; i<n; i++)
301 {
302 const MSignalPix &pix = (*this)[i];
303 if (!pix.IsPixelUsed())
304 continue;
305
306 Float_t testval = pix.GetErrorPhot();
307
308 if (geom)
309 testval *= geom->GetPixRatio(i);
310
311 if (testval < minval)
312 minval = testval;
313 }
314 return minval;
315}
316
317// --------------------------------------------------------------------------
318//
319// get the maximum ratio of photons/error
320// If you specify a geometry the number of photons is weighted with the
321// area of the pixel
322//
323Float_t MSignalCam::GetErrorPhotMax(const MGeomCam *geom) const
324{
325 const UInt_t n = GetNumPixels();
326
327 if (n <= 0)
328 return 50.;
329
330 Float_t maxval = -FLT_MAX;
331
332 for (UInt_t i=0; i<n; i++)
333 {
334 const MSignalPix &pix = (*this)[i];
335 if (!pix.IsPixelUsed())
336 continue;
337
338 Float_t testval = pix.GetErrorPhot();
339
340 if (geom)
341 testval *= geom->GetPixRatio(i);
342
343 if (testval > maxval)
344 maxval = testval;
345 }
346 return maxval;
347}
348
349MSignalPix *MSignalCam::AddPixel(Int_t idx, Float_t nph, Float_t er)
350{
351 MSignalPix *pix = static_cast<MSignalPix*>((*fPixels)[idx]);
352 pix->Set(nph, er);
353 return pix;
354}
355
356// --------------------------------------------------------------------------
357//
358// This function recursively finds all pixels of one island and assigns
359// the number num as island number to the pixel.
360//
361// 1) Check whether a pixel with the index idx exists, is unused
362// and has not yet a island number assigned.
363// 2) Assign the island number num to the pixel
364// 3) Loop over all its neighbors taken from the geometry geom. For all
365// neighbors recursively call this function (CalcIsland)
366// 4) Sum the size of the pixel and all neighbors newly assigned
367// (by CalcIsland) to this island
368//
369// Returns the sum of the pixel size.
370//
371Double_t MSignalCam::CalcIsland(const MGeomCam &geom, Int_t idx, Int_t num)
372{
373 // Get the pixel information of a pixel with this index
374 MSignalPix &pix = (*this)[idx];
375
376 // If an island number was already assigned to this pixel... do nothing.
377 if (pix.GetIdxIsland()>=0)
378 return 0;
379
380 // If the pixel is an unused pixel... do nothing.
381 if (!pix.IsPixelUsed())
382 return 0;
383
384 // Assign the new island number num to this used pixel
385 pix.SetIdxIsland(num);
386
387 // Get the geometry information (neighbors) of this pixel
388 const MGeom &gpix = geom[idx];
389
390 // Get the size of this pixel
391 Double_t size = pix.GetNumPhotons();
392
393 // Now do the same with all its neighbors and sum the
394 // sizes which they correspond to
395 const Int_t n = gpix.GetNumNeighbors();
396 for (int i=0; i<n; i++)
397 size += CalcIsland(geom, gpix.GetNeighbor(i), num);
398
399 // return size of this (sub)cluster
400 return size;
401}
402
403// --------------------------------------------------------------------------
404//
405// Each pixel which is maked as used is assigned an island number
406// (starting from 0). A pixel without an island number assigned
407// has island number -1.
408//
409// The index 0 corresponds to the island with the highest size (sum
410// of GetNumPhotons() in island). The size is decreasing with
411// increasing indices.
412//
413// The information about pixel neighbory is taken from the geometry
414// MGeomCam geom;
415//
416// You can access this island number of a pixel with a call to
417// MSignalPix->GetIdxIsland. The total number of islands available
418// can be accessed with MSignalCam->GetNumIslands.
419//
420// CalcIslands returns the number of islands found. If an error occurs,
421// eg the geometry has less pixels than the highest index stored, -1 is
422// returned.
423//
424Int_t MSignalCam::CalcIslands(const MGeomCam &geom)
425{
426 const UInt_t numpix = GetNumPixels();
427
428 if (/*fMaxIndex<0 ||*/ numpix==0)
429 {
430 *fLog << err << "ERROR - MSignalCam doesn't contain pixels!" << endl;
431 fNumIslands = 0;
432 return -1;
433 }
434/*
435 if ((UInt_t)fMaxIndex>=geom.GetNumPixels())
436 {
437 *fLog << err << "ERROR - MSignalCam::CalcIslands: Size mismatch - geometry too small!" << endl;
438 return -1;
439 }
440 */
441 // Create a list to hold the sizes of the islands (The maximum
442 // number of islands possible is roughly fNumPixels/4)
443 TArrayD size(numpix/3);
444
445 // Calculate Islands
446 Int_t n=0;
447 Float_t totsize = 0;
448
449 for (UInt_t idx=0; idx<numpix; idx++)
450 {
451 const MSignalPix &pix = (*this)[idx];
452 if (!pix.IsPixelUsed())
453 continue;
454
455 // This 'if' is necessary (although it is done in GetIsland, too)
456 // because otherwise the counter (n) would be wrong.
457 // So only 'start' a new island for used pixels (selected by
458 // using the Iterator) which do not yet belong to another island.
459 if (pix.GetIdxIsland()<0)
460 {
461 // Don't put this in one line! n is used twice...
462 const Double_t sz = CalcIsland(geom, idx/*pix->GetPixId()*/, n);
463 size[n++] = sz;
464 totsize += sz;
465 }
466 }
467
468 // Create an array holding the indices
469 TArrayI idxarr(n);
470
471 // Sort the sizes descending
472 TMath::Sort(n, size.GetArray(), idxarr.GetArray(), kTRUE);
473
474 // Replace island numbers by size indices -- After this
475 // islands indices are sorted by the island size
476 for (UInt_t idx=0; idx<numpix; idx++)
477 {
478 MSignalPix &pix = (*this)[idx];
479 if (!pix.IsPixelUsed())
480 continue;
481
482 const Short_t i = pix.GetIdxIsland();
483
484 // Find new index
485 Short_t j;
486 for (j=0; j<n; j++)
487 if (idxarr[j]==i)
488 break;
489
490 pix.SetIdxIsland(j==n ? -1 : j);
491 }
492
493 // Now assign number of islands found
494 fNumIslands = n;
495 fSizeSubIslands = n>0 ? totsize-size[idxarr[0]] : 0;
496 fSizeMainIsland = n>0 ? size[idxarr[0]] : 0;
497
498 // return number of island
499 return fNumIslands;
500}
501
502// --------------------------------------------------------------------------
503//
504// Compares the cleaning (fRing and fIsCore) of this object and the
505// argument. Return the result (kFALSE if different).
506//
507// If differences are found the pixels (the pixel from this object first)
508// is printed.
509//
510Bool_t MSignalCam::CompareCleaning(const MSignalCam &cam) const
511{
512 const UInt_t n = GetNumPixels();
513
514 if (n != cam.GetNumPixels())
515 {
516 *fLog << warn << "MSignalCam::CompareCleaning - Number of pixels mismatch." << endl;
517 return kFALSE;
518 }
519
520 for (UInt_t i=0; i<n; i++)
521 {
522 const MSignalPix &p = (*this)[i];
523 const MSignalPix &q = cam[i];
524
525 if (p.GetRing()==q.GetRing() && p.IsPixelCore()==q.IsPixelCore())
526 continue;
527
528 *fLog << warn << "MSignalCam::CompareCleaning - Pixel #" << i << " mismatch." << endl;
529 p.Print();
530 q.Print();
531 return kFALSE;
532 }
533 return kTRUE;
534}
535
536// --------------------------------------------------------------------------
537//
538// Compares the islands (fIdxIsland) of this object and the
539// argument. Return the result (kFALSE if different).
540//
541// If differences are found the pixels (the pixel from this object first)
542// is printed.
543//
544Bool_t MSignalCam::CompareIslands(const MSignalCam &cam) const
545{
546 const UInt_t n = GetNumPixels();
547
548 if (n != cam.GetNumPixels())
549 {
550 *fLog << warn << "MSignalCam::CompareIslands - Number of pixels mismatch." << endl;
551 return kFALSE;
552 }
553
554 for (UInt_t i=0; i<n; i++)
555 {
556 const MSignalPix &p = (*this)[i];
557 const MSignalPix &q = cam[i];
558
559 if (p.GetIdxIsland()==q.GetIdxIsland())
560 continue;
561
562 *fLog << warn << "MSignalCam::CompareIslands - Pixel #" << i << " mismatch." << endl;
563 p.Print();
564 q.Print();
565 return kFALSE;
566 }
567
568 return kTRUE;
569}
570
571// --------------------------------------------------------------------------
572//
573// Returns, depending on the type flag:
574//
575// 0: Number of Photons*PixRatio <default>
576// 1: Error*sqrt(PixRatio)
577// 2: Cleaning level = Num Photons*sqrt(PixRatio)/Error
578// 3: Number of Photons
579// 4: Error
580// 5: Island index
581// 6: arrival time of mapped pixels
582// 7: arrival time if signa avove 20phe
583// 8: arrival time
584// 10: as 0, but returns kFALSE if signal <=0
585// 11: as 8, but returns kFALSE if signal <=0
586//
587Bool_t MSignalCam::GetPixelContent(Double_t &val, Int_t idx, const MGeomCam &cam, Int_t type) const
588{
589 if (idx<0 || (UInt_t)idx>GetNumPixels())
590 return kFALSE;
591
592 const MSignalPix &pix = (*this)[idx];
593
594 // Used inlcudes status unampped
595 if (!pix.IsPixelUsed() && (type<6 || type==8))
596 return kFALSE;
597
598 const Double_t ratio = cam.GetPixRatio(idx);
599
600 switch (type)
601 {
602 case 1: // scaled error of phtoons
603 val = pix.GetErrorPhot()*TMath::Sqrt(ratio);
604 return kTRUE;
605
606 case 2: // significance of number of photons
607 if (pix.GetErrorPhot()<=0)
608 return kFALSE;
609 val = pix.GetNumPhotons()*TMath::Sqrt(ratio)/pix.GetErrorPhot();
610 return kTRUE;
611
612 case 3: // number of photo electrons
613 val = pix.GetNumPhotons();
614 break;
615
616 case 4: // error of signal
617 val = pix.GetErrorPhot();
618 break;
619
620 case 5: // index of island
621 val = pix.GetIdxIsland();
622 break;
623
624 case 6: // arrival time of mapped pixels only
625 if (pix.IsPixelUnmapped())
626 return kFALSE;
627 val = pix.GetArrivalTime();
628 break;
629
630 case 7: // pulse position above 50phe
631 // The number of photons is not scaled with the ratio because
632 // otherwise to many large pixels survive (maybe because the
633 // fluctuations scale different than expected)
634 if (pix.IsPixelUnmapped() || pix.GetNumPhotons()<50)
635 return kFALSE;
636 val = pix.GetArrivalTime();
637 break;
638
639 case 8: // arrival time
640 val = pix.GetArrivalTime();
641 break;
642
643 /*
644 case 10: // lo gain time
645 if (pix.IsPixelUnmapped() || !pix.IsLoGainUsed() ||
646 pix.GetNumPhotons()<320)
647 return kFALSE;
648 val = pix.GetArrivalTime();
649 return kTRUE;
650
651 case 11: // hi gain time
652 // The number of photons is not scaled with the ratio because
653 // otherwise to many large pixels survive (maybe because the
654 // fluctuations scale different than expected)
655 if (pix.IsPixelUnmapped() || pix.IsLoGainUsed() ||
656 pix.GetNumPhotons()<50)
657 return kFALSE;
658 val = pix.GetArrivalTime();
659 return kTRUE;
660 */
661
662 case 10:
663 val = pix.GetNumPhotons()*ratio;
664 return val>0;
665
666 case 11:
667 val = pix.GetArrivalTime();
668 return pix.GetNumPhotons()>0;
669
670 case 9:
671 default:
672 val = pix.GetNumPhotons()*ratio;
673 return kTRUE;
674 }
675 return kTRUE;
676}
677
678void MSignalCam::DrawPixelContent(Int_t num) const
679{
680 *fLog << warn << "MSignalCam::DrawPixelContent - not available." << endl;
681}
682
683TObject *MSignalCamIter::Next()
684{
685 if (!fUsedOnly)
686 {
687 fIdx++;
688 return TObjArrayIter::Next();
689 }
690
691 MSignalPix *pix;
692 while ((pix = (MSignalPix*)TObjArrayIter::Next()))
693 {
694 fIdx++;
695 if (pix->IsPixelUsed())
696 return pix;
697 }
698 return pix;
699}
Note: See TracBrowser for help on using the repository browser.