source: branches/Mars_McMismatchStudy/mimage/MImgCleanTime.cc@ 17985

Last change on this file since 17985 was 17860, checked in by tbretz, 11 years ago
Fixed a compilation problem.
File size: 7.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, 4/2014 <mailto:tbretz@phys.ethz.ch>
19!
20! Copyright: MAGIC Software Development, 2014
21!
22!
23\* ======================================================================== */
24
25/////////////////////////////////////////////////////////////////////////////
26//
27// MImgCleanTime
28//
29// Input Containers:
30// MGeomCam
31// MSignalCam
32//
33// Output Containers:
34// MSignalCam
35//
36/////////////////////////////////////////////////////////////////////////////
37#include "MImgCleanTime.h"
38
39#include <TEnv.h>
40
41#include "MLog.h"
42#include "MLogManip.h"
43
44#include "MParList.h"
45
46#include "MGeomPix.h"
47#include "MGeomCam.h"
48
49#include "MSignalPix.h"
50#include "MSignalCam.h"
51
52ClassImp(MImgCleanTime);
53
54using namespace std;
55
56// --------------------------------------------------------------------------
57//
58// Default constructor. Here you can specify the cleaning method and levels.
59// If you don't specify them the 'common standard' values 3.0 and 2.5 (sigma
60// above mean) are used.
61// Here you can also specify how many rings around the core pixels you want
62// to analyze (with the fixed lvl2). The default value for "rings" is 1.
63//
64MImgCleanTime::MImgCleanTime(const char *name, const char *title)
65 : fCam(0), fEvt(0), fMinCount(0), fMinSize(25), fDeltaT(2.5*17.5*0.1111),
66 fNameSignalCam("MSignalCam")
67{
68 fName = name ? name : "MImgCleanTime";
69 fTitle = title ? title : "Task to perform image cleaning";
70}
71
72// --------------------------------------------------------------------------
73//
74// Check if MEvtHeader exists in the Parameter list already.
75// if not create one and add them to the list
76//
77Int_t MImgCleanTime::PreProcess (MParList *pList)
78{
79 fCam = (MGeomCam*)pList->FindObject("MGeomCam");
80 if (!fCam)
81 {
82 *fLog << err << "MGeomCam not found (no geometry information available)... aborting." << endl;
83 return kFALSE;
84 }
85 fEvt = (MSignalCam*)pList->FindObject(fNameSignalCam, "MSignalCam");
86 if (!fEvt)
87 {
88 *fLog << err << fNameSignalCam << " [MSignalCam] not found... aborting." << endl;
89 return kFALSE;
90 }
91
92 return kTRUE;
93}
94
95Island MImgCleanTime::CalcIsland(MSignalPix &pix1, const MGeom &gpix1, const uint16_t &island1)
96{
97 pix1.SetIdxIsland(island1);
98
99 const Float_t time1 = pix1.GetArrivalTime();
100
101 Island island(pix1.GetNumPhotons(), time1);
102
103 for (UInt_t i=0; i<gpix1.GetNumNeighbors(); i++)
104 {
105 const Int_t idx2 = gpix1.GetNeighbor(i);
106
107 MSignalPix &pix2 = (*fEvt)[idx2];
108
109 const Int_t island2 = pix2.GetIdxIsland();
110 if (island2>=0)
111 {
112 if (island1==island2)
113 continue;
114
115 // The entries are sorted naturally, therefore
116 // there is no need to check the whole array
117 const auto it = fContacts.rbegin();
118 if (it==fContacts.rend() ||
119 (it->first!=island1 && it->second!=island2))
120 fContacts.emplace_back(island1, island2);
121
122 continue;
123 }
124
125 const Float_t time2 = pix2.GetArrivalTime();
126 if (fabs(time2-time1)<fDeltaT) // FIXME: Scale with distance?
127 island += CalcIsland(pix2, (*fCam)[idx2], island1);
128 }
129
130 return island;
131}
132
133// --------------------------------------------------------------------------
134//
135// Cleans the image.
136//
137Int_t MImgCleanTime::Process()
138{
139 // The assumption is that all pixels contain valid data
140 const UInt_t npix = fEvt->GetNumPixels();
141 for (UInt_t i=0; i<npix; i++)
142 {
143 MSignalPix &pix = (*fEvt)[i];
144 pix.SetPixelUnused();
145 pix.SetPixelCore(kFALSE);
146 pix.SetIdxIsland(-1);
147 }
148
149 // Contains the id of the island for each pixel
150 fIslands.clear();
151 fContacts.clear();
152
153 // Start with island idx==0
154 UShort_t idx = 0;
155 for (UInt_t i=0; i<npix; i++)
156 {
157 MSignalPix &spix = (*fEvt)[i];
158
159 // The following might be much more efficient and faster
160 // if we omit small counted and sized islands already,
161 // but what is the disadvantage?
162 if (spix.GetIdxIsland()<0)
163 fIslands.emplace_back(CalcIsland(spix, (*fCam)[i], idx++));
164 }
165
166 fLut.resize(fIslands.size());
167 for (UInt_t i=0; i<fLut.size(); i++)
168 fLut[i] = i;
169
170 // Unify touching islands if their arrival time ranges overlap
171 for (auto it=fContacts.begin(); it!=fContacts.end(); it++)
172 {
173 const uint16_t &idx1 = fLut[it->first];
174 const uint16_t &idx2 = fLut[it->second];
175 if (idx1==idx2)
176 continue;
177
178 Island &I1 = fIslands[idx1];
179 const Island &I2 = fIslands[idx2];
180
181 // FIXME: Put a limit on count? Put a limit on size?
182 // The ideal cut would be to allow that a single
183 // pixel still can connect two clusters
184 if (I1.count==1 || I1.count==1)
185 continue;
186
187 if (I1.max<I2.min || I2.max<I1.min)
188 continue;
189
190 // Combine both islands
191 I1 += I2;
192
193 // Remove idx2 from the list
194 for (auto is=fLut.begin(); is!=fLut.end(); is++)
195 if (*is==idx2)
196 *is=idx1;
197
198 // Erase entry form contact list
199 fContacts.erase(it);
200
201 // Start over again
202 it = fContacts.begin();
203 }
204
205 uint16_t num_islands = 0;
206 double size_main = 0;
207 double size_tot = 0;
208
209 vector<bool> used(fIslands.size());
210 for (UInt_t i=0; i<npix; i++)
211 {
212 MSignalPix &pix = (*fEvt)[i];
213
214 // At the end every pixel has an island assigned
215 const Short_t ii = fLut[pix.GetIdxIsland()];
216
217 const Island &island = fIslands[ii];
218
219 if (island.count<=fMinCount || island.size<=fMinSize)
220 continue;
221
222 if (!used[ii])
223 {
224 used[ii] = true;
225 num_islands++;
226
227 size_tot += island.size;
228 if (island.size>size_main)
229 size_main = island.size;
230 }
231
232 pix.SetPixelUsed();
233 pix.SetPixelCore();
234 }
235
236 fEvt->SetIslandInfo(num_islands, size_main, size_tot-size_main);
237
238 return kTRUE;
239}
240
241// --------------------------------------------------------------------------
242//
243Int_t MImgCleanTime::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
244{
245 Bool_t rc = kFALSE;
246 if (IsEnvDefined(env, prefix, "MinSize", print))
247 {
248 rc = kTRUE;
249 SetMinSize(GetEnvValue(env, prefix, "MinSize", fMinSize));
250 }
251 if (IsEnvDefined(env, prefix, "MinCount", print))
252 {
253 rc = kTRUE;
254 SetMinSize(GetEnvValue(env, prefix, "MinCount", (Int_t)fMinCount));
255 }
256 if (IsEnvDefined(env, prefix, "DeltaT", print))
257 {
258 rc = kTRUE;
259 SetMinSize(GetEnvValue(env, prefix, "DeltaT", fDeltaT));
260 }
261
262 return rc;
263}
Note: See TracBrowser for help on using the repository browser.