source: trunk/Mars/msignal/MSignalCam.cc @ 18243

Last change on this file since 18243 was 18243, checked in by dorner, 5 years ago
added new parameter fTimeSlope to the output of the calibration
File size: 18.4 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// 12: time slope
587//
588Bool_t MSignalCam::GetPixelContent(Double_t &val, Int_t idx, const MGeomCam &cam, Int_t type) const
589{
590    if (idx<0 || (UInt_t)idx>GetNumPixels())
591        return kFALSE;
592
593    const MSignalPix &pix = (*this)[idx];
594
595    // Used inlcudes status unampped
596    if (!pix.IsPixelUsed() && (type<6 || type==8))
597        return kFALSE;
598
599    const Double_t ratio = cam.GetPixRatio(idx);
600
601    switch (type)
602    {
603    case 1: // scaled error of phtoons
604        val = pix.GetErrorPhot()*TMath::Sqrt(ratio);
605        return kTRUE;
606
607    case 2: // significance of number of photons
608        if (pix.GetErrorPhot()<=0)
609            return kFALSE;
610        val = pix.GetNumPhotons()*TMath::Sqrt(ratio)/pix.GetErrorPhot();
611        return kTRUE;
612
613    case 3: // number of photo electrons
614        val = pix.GetNumPhotons();
615        break;
616
617    case 4: // error of signal
618        val = pix.GetErrorPhot();
619        break;
620
621    case 5: // index of island
622        val = pix.GetIdxIsland();
623        break;
624
625    case 6: // arrival time of mapped pixels only
626        if (pix.IsPixelUnmapped())
627            return kFALSE;
628        val = pix.GetArrivalTime();
629        break;
630
631    case 7: // pulse position above 50phe
632        // The number of photons is not scaled with the ratio because
633        // otherwise to many large pixels survive (maybe because the
634        // fluctuations scale different than expected)
635        if (pix.IsPixelUnmapped() || pix.GetNumPhotons()<50)
636            return kFALSE;
637        val = pix.GetArrivalTime();
638        break;
639
640    case 8: // arrival time
641        val = pix.GetArrivalTime();
642        break;
643
644        /*
645    case 10: // lo gain time
646        if (pix.IsPixelUnmapped() || !pix.IsLoGainUsed() ||
647            pix.GetNumPhotons()<320)
648            return kFALSE;
649        val = pix.GetArrivalTime();
650        return kTRUE;
651
652    case 11: // hi gain time
653        // The number of photons is not scaled with the ratio because
654        // otherwise to many large pixels survive (maybe because the
655        // fluctuations scale different than expected)
656        if (pix.IsPixelUnmapped() || pix.IsLoGainUsed() ||
657            pix.GetNumPhotons()<50)
658            return kFALSE;
659        val = pix.GetArrivalTime();
660        return kTRUE;
661        */
662
663    case 10:
664        val = pix.GetNumPhotons()*ratio;
665        return val>0;
666
667    case 11:
668        val = pix.GetArrivalTime();
669        return pix.GetNumPhotons()>0;
670
671    case 12: // time slope
672        val = pix.GetTimeSlope();
673        break;
674
675    case 9:
676    default:
677        val = pix.GetNumPhotons()*ratio;
678        return kTRUE;
679    }
680    return kTRUE;
681}
682
683void MSignalCam::DrawPixelContent(Int_t num) const
684{
685    *fLog << warn << "MSignalCam::DrawPixelContent - not available." << endl;
686}
687
688TObject *MSignalCamIter::Next()
689{
690    if (!fUsedOnly)
691    {
692        fIdx++;
693        return TObjArrayIter::Next();
694    }
695
696    MSignalPix *pix;
697    while ((pix = (MSignalPix*)TObjArrayIter::Next()))
698    {
699        fIdx++;
700        if (pix->IsPixelUsed())
701            return pix;
702    }
703    return pix;
704}
Note: See TracBrowser for help on using the repository browser.