source: branches/MarsMoreSimulationTruth/mtemp/mifae/library/MIslandsClean.cc

Last change on this file was 5394, checked in by aliu, 20 years ago
*** empty log message ***
File size: 9.0 KB
Line 
1/* ======================================================================== *\
2!
3!
4! Author(s): Ester Aliu, 3/2004
5!
6!
7\* ======================================================================== */
8
9/////////////////////////////////////////////////////////////////////////////
10//
11// MIslandClean
12//
13// The Island Cleaning task selects the islands you use for the Hillas
14// parameters calculation, after the normal image cleaning.
15//
16// There are two methods to make the selection:
17//
18// - No time method, as used is Whipple. It calculates an island parameter
19// called "signal to noise" and adds a new threshold that eliminates some
20// of the islands. The way the island is eliminated is seeting the pixels
21// of the islands as UNUSED
22//
23// - Time method, taking profit of the time information in MAGIC.
24// Calculates the maximum time difference (in time slices) for each island
25// and corrects for the island size. With an other new threshold "noise"
26// islands are supposed to be eliminated.
27//
28// Other cleanings that are allowed in this code are:
29//
30// - Resting only with the continent, i.e. the larger island
31//
32// Example:
33//
34// MIslands isl;
35// isl.SetName("MIslands1");
36
37// MImgCleanStd clean;
38// MIslandClean islclean(0.2);
39// islclean.SetInputName("MIslands1");
40// islclean.SetMethod(0); // for timing method
41//
42// tlist.AddToList(&clean);
43// tlist.AddToList(&islclean);
44//
45//
46// Input Containers:
47// MGeomCam
48// MCerPhotEvt
49// MPedestalCam
50// MArrivalTime
51// MIslands
52//
53// Output Containers:
54// MCerPhotEvt
55//
56/////////////////////////////////////////////////////////////////////////////
57#include "MIslandsClean.h"
58
59#include <stdlib.h> // atof
60#include <fstream> // ofstream, SavePrimitive
61
62#include "MLog.h"
63#include "MLogManip.h"
64
65#include "MIslands.h"
66#include "MImgIsland.h"
67
68#include "MParList.h"
69
70#include "MGeomPix.h"
71#include "MGeomCam.h"
72
73#include "MCerPhotPix.h"
74#include "MCerPhotEvt.h"
75
76#include "MPedestalCam.h"
77#include "MPedestalPix.h"
78
79#include "MArrivalTimeCam.h"
80#include "MArrivalTimePix.h"
81
82ClassImp(MIslandsClean);
83
84
85using namespace std;
86
87// --------------------------------------------------------------------------
88//
89// Default constructor.
90//
91MIslandsClean::MIslandsClean(const Float_t newThres, const char *name, const char *title)
92 : fIsl(NULL), fIslandCleaningMethod(kNoTiming), fIslCleanThres(newThres)
93{
94 fName = name ? name : "MIslandsClean";
95 fTitle = title ? title : "Clean islands";
96}
97
98
99Int_t MIslandsClean::PreProcess (MParList *pList)
100{
101 fCam = (MGeomCam*)pList->FindObject(AddSerialNumber("MGeomCam"));
102 if (!fCam)
103 {
104 *fLog << dbginf << "MGeomCam not found (no geometry information available)... aborting." << endl;
105 return kFALSE;
106 }
107
108 fEvt = (MCerPhotEvt*)pList->FindObject(AddSerialNumber("MCerPhotEvt"));
109 if (!fEvt)
110 {
111 *fLog << dbginf << "MCerPhotEvt not found... aborting." << endl;
112 return kFALSE;
113 }
114
115 fPed = (MPedestalCam*)pList->FindObject(AddSerialNumber("MPedestalCam"));
116 if (!fPed)
117 {
118 *fLog << dbginf << "MPedestalCam not found... aborting." << endl;
119 return kFALSE;
120 }
121
122 fTime = (MArrivalTimeCam*)pList->FindObject(AddSerialNumber("MArrivalTimeCam"));
123 if (!fTime)
124 {
125 *fLog << dbginf << "MArrivalTimeCam not found... aborting." << endl;
126 return kFALSE;
127 }
128
129 if (strlen(fIslName) > 0)
130 fIsl = (MIslands*)pList->FindObject(AddSerialNumber(fIslName));
131 else
132 fIsl = (MIslands*)pList->FindObject(AddSerialNumber("MIslands"));
133 if (!fIsl)
134 {
135 *fLog << dbginf << "MIslands not found... aborting." << endl;
136 return kFALSE;
137 }
138
139 return kTRUE;
140}
141
142
143
144Int_t MIslandsClean::Process()
145{
146
147 MImgIsland *imgIsl = new MImgIsland;
148 TIter Next(fIsl->GetList());
149
150 Int_t pixNum = 0;
151 Int_t idPix = -1;
152
153 ////////////////////////////////////////////////////
154 //
155 // TIME SPREAD ISLAND CLEANING
156 // eliminates the island with a time spread
157 // higher than a given limit
158 ///////////////////////////////////////////////////
159 if( fIslandCleaningMethod == 0 ){
160 while ((imgIsl=(MImgIsland*)Next())) {
161
162 //fIslCleanThreshold has different values, FIXME, put two variables
163 if(imgIsl->GetTimeSpread() > fIslCleanThres)
164 {
165 pixNum = imgIsl->GetPixNum();
166
167 for(Int_t k = 0; k<pixNum; k++)
168 {
169 idPix = imgIsl->GetPixList(k);
170 MCerPhotPix &pix = (*fEvt)[idPix];
171 pix.SetPixelUnused();
172 }
173 }
174 }
175
176 }
177
178
179 ////////////////////////////////////////////////////
180 //
181 // SIGNAL TO NOISE ISLAND CLEANING
182 // eliminates the island with a signal-to-noise
183 // lower than a given limit
184 ///////////////////////////////////////////////////
185 else if ( fIslandCleaningMethod == 1 ) {
186 while ((imgIsl=(MImgIsland*)Next())) {
187
188 if(imgIsl->GetSigToNoise() < fIslCleanThres)
189 {
190 pixNum = imgIsl->GetPixNum();
191
192 for(Int_t k = 0; k<pixNum; k++)
193 {
194 idPix = imgIsl->GetPixList(k);
195 MCerPhotPix &pix = (*fEvt)[idPix];
196 pix.SetPixelUnused();
197 }
198 }
199 }
200 }
201
202
203 ////////////////////////////////////////////////////
204 //
205 // ISLAND SIZE OVER EVENT SIZE ISLAND CLEANING
206 // eliminates the island with an island size over event size
207 // lower than a given limit
208 ///////////////////////////////////////////////////
209 else if ( fIslandCleaningMethod == 2 ) {
210 Float_t size = 0;
211 while ((imgIsl=(MImgIsland*)Next())) {
212 size += imgIsl->GetSizeIsl();
213 }
214
215 Next.Reset();
216 while ((imgIsl=(MImgIsland*)Next())) {
217
218 if(imgIsl->GetSizeIsl()/size < fIslCleanThres)
219 {
220 pixNum = imgIsl->GetPixNum();
221
222 for(Int_t k = 0; k<pixNum; k++)
223 {
224 idPix = imgIsl->GetPixList(k);
225 MCerPhotPix &pix = (*fEvt)[idPix];
226 pix.SetPixelUnused();
227 }
228 }
229 }
230 }
231
232
233 ////////////////////////////////////////////////////
234 //
235 // CONTINENT ISLAND CLEANING
236 // eliminates all the islands except the continent
237 // i.e. the larger island in the event
238 // according the number of pixels
239 ///////////////////////////////////////////////////
240 else if( fIslandCleaningMethod == 3 ){
241 Int_t i = 0;
242 while ((imgIsl=(MImgIsland*)Next())) {
243 if (i != 0){
244
245 pixNum = imgIsl->GetPixNum();
246
247 for(Int_t k = 0; k<pixNum; k++)
248 {
249 idPix = imgIsl->GetPixList(k);
250 MCerPhotPix &pix = (*fEvt)[idPix];
251 pix.SetPixelUnused();
252 }
253 }
254 i++;
255 }
256 }
257
258
259 ////////////////////////////////////////////////////
260 //
261 // LARGER and SECOND LARGER ISLAND CLEANING
262 // eliminates all the islands except the two biggest
263 // ones according size
264 //
265 ///////////////////////////////////////////////////
266 else if( fIslandCleaningMethod == 4 ){
267
268 Int_t i = 0;
269 Int_t islnum = fIsl->GetIslNum();
270 Float_t islSize[islnum];
271 Int_t islIdx[islnum];
272
273 for (Int_t j = 0; j<islnum ; j++){
274 islSize[j] = -1;
275 islIdx[j] = -1;
276 }
277
278 while ((imgIsl=(MImgIsland*)Next())) {
279
280 islSize[i] = imgIsl->GetSizeIsl();
281 i++;
282 }
283
284
285 TMath::Sort(islnum, islSize, islIdx, kTRUE);
286
287 i = 0;
288 Next.Reset();
289 while ((imgIsl=(MImgIsland*)Next())) {
290 if (islnum > 1 && islIdx[0]!=i && islIdx[1]!=i){
291
292 //cout << "removed " << i << " isl 0" << islIdx[0] << " isl 1" << islIdx[1] << endl;
293
294 pixNum = imgIsl->GetPixNum();
295
296 for(Int_t k = 0; k<pixNum; k++)
297 {
298 idPix = imgIsl->GetPixList(k);
299 MCerPhotPix &pix = (*fEvt)[idPix];
300 pix.SetPixelUnused();
301 }
302 }
303 i++;
304 }
305 }
306
307
308
309 ///////////////////////////////////////////////////////
310 //
311 // LARGER and SECOND LARGER ISLAND CLEANING II
312 // eliminates all the islands except the two biggest
313 // ones according size, if the second one almost has
314 // the 80% of the size of the biggest one
315 //
316 //
317 //////////////////////////////////////////////////////
318 else if( fIslandCleaningMethod == 5 ){
319
320 Int_t i = 0;
321 Int_t islnum = fIsl->GetIslNum();
322 Float_t islSize[islnum];
323 Int_t islIdx[islnum];
324
325 for (Int_t j = 0; j<islnum ; j++){
326 islSize[j] = -1;
327 islIdx[j] = -1;
328 }
329
330 while ((imgIsl=(MImgIsland*)Next())) {
331
332 islSize[i] = imgIsl->GetSizeIsl();
333 i++;
334 }
335
336 TMath::Sort(islnum, islSize, islIdx, kTRUE);
337
338 i = 0;
339 Next.Reset();
340
341 while ((imgIsl=(MImgIsland*)Next())) {
342
343 if (i!=0 && islIdx[0]!=i && islIdx[1]!=i){
344
345 pixNum = imgIsl->GetPixNum();
346
347 for(Int_t k = 0; k<pixNum; k++)
348 {
349 idPix = imgIsl->GetPixList(k);
350 MCerPhotPix &pix = (*fEvt)[idPix];
351 pix.SetPixelUnused();
352 }
353 }
354 else if(i!=0 && islSize[islIdx[i]]<0.6*islSize[islIdx[0]]){
355
356 pixNum = imgIsl->GetPixNum();
357
358 for(Int_t k = 0; k<pixNum; k++)
359 {
360 idPix = imgIsl->GetPixList(k);
361 MCerPhotPix &pix = (*fEvt)[idPix];
362 pix.SetPixelUnused();
363 }
364 }
365 i++;
366 }
367 }
368
369 fEvt->SetReadyToSave();
370
371 return kTRUE;
372
373}
Note: See TracBrowser for help on using the repository browser.