source: trunk/MagicSoft/Mars/manalysis/MBlindPixelCalc.cc@ 2908

Last change on this file since 2908 was 2502, checked in by tbretz, 21 years ago
*** empty log message ***
File size: 12.9 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): Oscar Blanch 12/2001 <mailto:blanch@ifae.es>
19! Author(s): Thomas Bretz 08/2002 <mailto:tbretz@astro.uni.wuerzburg.de>
20!
21! Copyright: MAGIC Software Development, 2000-2003
22!
23!
24\* ======================================================================== */
25
26/////////////////////////////////////////////////////////////////////////////
27//
28// MBlindPixelCalc
29//
30// This is the specific image cleaning for a list of pixels. This task
31// sets the pixels listed in fPixelsIdx to unused so they should not be
32// used for analysis (eg calculation of hillas parameters).
33//
34// You can use MBlindPixelCalc::SetUseInterpolation to replaced the
35// blind pixels by the average of its neighbors instead of unmapping
36// them. If you want to include the central pixel use
37// MBlindPixelCalc::SetUseCentralPixel.
38//
39// You have three options:
40// 1) Call SetUseBlindPixels():
41// This will take an existing MBlindPixels container filled from
42// elsewhere (eg. MCT1ReadPreProc) and use this pixels as blind
43// pixels.
44// 2) Call SetPixels():
45// This will setup an array with pixel numbers. These pixels are used
46// as blind pixels.
47// 3) Neither 1) nor 2)
48// This options tries to identify the starfield from the
49// MMcRunHeader container and tries to identifies it. If it is known
50// (eg. Crab) the fixed build in pixel numbers are used as blind
51// pixels.
52//
53// If neither an array of pixels is given (or its size is 0) and
54// MMcRunHeader couldn't be found the task removes itself from the
55// tasklist.
56//
57// Implemented star fields:
58// - Crab: 400, 401, 402, 437, 438, 439
59//
60// Input Containers:
61// MCerPhotEvt
62// [MBlindPixels]
63//
64// Output Containers:
65// MCerPhotEvt
66// MBlindPixels
67//
68/////////////////////////////////////////////////////////////////////////////
69#include "MBlindPixelCalc.h"
70
71#include <fstream>
72
73#include "MLog.h"
74#include "MLogManip.h"
75
76#include "MParList.h"
77
78#include "MGeomPix.h"
79#include "MGeomCam.h"
80#include "MCerPhotPix.h"
81#include "MCerPhotEvt.h"
82#include "MPedestalPix.h"
83#include "MPedestalCam.h"
84#include "MBlindPixels.h"
85
86#include "MMcRunHeader.hxx"
87
88ClassImp(MBlindPixelCalc);
89
90using namespace std;
91
92static const TString gsDefName = "MBlindPixelCalc";
93static const TString gsDefTitle = "Find hot spots (star, broken pixels, etc)";
94
95// --------------------------------------------------------------------------
96//
97// Default constructor.
98//
99MBlindPixelCalc::MBlindPixelCalc(const char *name, const char *title)
100 : fFlags(0)
101{
102 fName = name ? name : gsDefName.Data();
103 fTitle = title ? title : gsDefTitle.Data();
104}
105
106// --------------------------------------------------------------------------
107//
108// - Try to find or create MBlindPixels in parameter list.
109// - get the MCerPhotEvt from the parlist (abort if missing)
110// - if no pixels are given by the user try to determin the starfield
111// from the monte carlo run header.
112//
113Int_t MBlindPixelCalc::PreProcess (MParList *pList)
114{
115 if (TESTBIT(fFlags, kUseBlindPixels))
116 fPixels = (MBlindPixels*)pList->FindObject(AddSerialNumber("MBlindPixels"));
117 else
118 fPixels = (MBlindPixels*)pList->FindCreateObj(AddSerialNumber("MBlindPixels"));
119 if (!fPixels)
120 return kFALSE;
121
122 fEvt = (MCerPhotEvt*)pList->FindObject(AddSerialNumber("MCerPhotEvt"));
123 if (!fEvt)
124 {
125 *fLog << err << dbginf << "MCerPhotEvt not found... aborting." << endl;
126 return kFALSE;
127 }
128 fPed = (MPedestalCam*)pList->FindObject(AddSerialNumber("MPedestalCam"));
129 if (!fEvt)
130 {
131 *fLog << err << dbginf << "MPedestalCam not found... aborting." << endl;
132 return kFALSE;
133 }
134
135 fGeomCam = (MGeomCam*)pList->FindObject(AddSerialNumber("MGeomCam"));
136 if (!fGeomCam)
137 *fLog << warn << dbginf << "No camera geometry available... can't use interpolation." << endl;
138
139 if (TESTBIT(fFlags, kUseBlindPixels))
140 return kTRUE;
141
142 const UShort_t size = fPixelsIdx.GetSize();
143
144 if (size == 0)
145 {
146 if (!pList->FindObject("MMcRunHeader"))
147 {
148 *fLog << warn << "Warning - Neither blind pixels are given nor a MMcRunHeader was found... removing MBlindPixelCalc from list." << endl;
149 return kSKIP;
150 }
151 return kTRUE;
152 }
153
154 // Set as blind pixels the global blind pixels, which are given
155 // through the macros
156
157 UShort_t numids = fPixelsIdx.GetSize();
158
159 for(Int_t i = 0; i<numids; i++)
160 fPixels->SetPixelBlind(fPixelsIdx[i]);
161
162 return kTRUE;
163}
164
165// --------------------------------------------------------------------------
166//
167// Replaces each pixel (signal, signal error, pedestal, pedestal rms)
168// by the average of its surrounding pixels.
169// If TESTBIT(fFlags, kUseCentralPixel) is set the central pixel is also
170// included.
171//
172// NT: Informations about the interpolated pedestals are added.
173// When the option Interpolated is called, the blind pixel with the new
174// values of signal and fluttuation is included in the calculation of
175// the Image Parameters.
176//
177void MBlindPixelCalc::Interpolate() const
178{
179 const UShort_t entries = fGeomCam->GetNumPixels();
180
181 //
182 // Create arrays
183 //
184 Double_t *nphot = new Double_t[entries];
185 Double_t *perr = new Double_t[entries];
186 Double_t *ped = new Double_t[entries];
187 Double_t *pedrms = new Double_t[entries];
188
189 //
190 // Loop over all pixels
191 //
192 for (UShort_t i=0; i<entries; i++)
193 {
194 //
195 // Check whether pixel with idx i is blind
196 //
197 if (!fPixels->IsBlind(i))
198 continue;
199
200 //
201 // Get a pointer to this pixel. If it is not yet existing
202 // create a new entry for this pixel in MCerPhotEvt
203 //
204 const MCerPhotPix *pix = fEvt->GetPixById(i);
205 if (!pix)
206 pix = fEvt->AddPixel(i, 0, 0);
207
208 //
209 // Get the corresponding geometry and pedestal
210 //
211 const MGeomPix &gpix = (*fGeomCam)[i];
212 const MPedestalPix &ppix = (*fPed)[i];
213
214 // Do Not-Use-Central-Pixel
215 const Bool_t nucp = !TESTBIT(fFlags, kUseCentralPixel);
216
217 Int_t num = nucp ? 0 : 1;
218
219 nphot[i] = nucp ? 0 : pix->GetNumPhotons();
220 ped[i] = nucp ? 0 : ppix.GetPedestal();
221 perr[i] = nucp ? 0 : Sqr(pix->GetErrorPhot());
222 pedrms[i] = nucp ? 0 : Sqr(ppix.GetPedestalRms());
223
224 //
225 // The values are rescaled to the small pixels area for the right comparison
226 //
227 const Double_t ratio = fGeomCam->GetPixRatio(i);
228
229 nphot[i] *= ratio;
230 ped[i] *= ratio;
231 perr[i] *= Sqr(ratio);
232 pedrms[i] *= Sqr(pedrms[i]);
233
234 //
235 // Loop over all its neighbors
236 //
237 const Int_t n = gpix.GetNumNeighbors();
238 for (int j=0; j<n; j++)
239 {
240 const UShort_t nidx = gpix.GetNeighbor(j);
241
242 //
243 // Do not use blind neighbors
244 //
245 if (fPixels->IsBlind(nidx))
246 continue;
247
248 //
249 // Check whether the neighbor has a signal stored
250 //
251 const MCerPhotPix *evtpix = fEvt->GetPixById(nidx);
252 if (!evtpix)
253 continue;
254
255 //
256 // Get the geometry for the neighbor
257 //
258 const Double_t nratio = fGeomCam->GetPixRatio(nidx);
259 MPedestalPix &nppix = (*fPed)[nidx];
260
261 //
262 //The error is calculated as the quadratic sum of the errors
263 //
264 ped[i] += nratio*nppix.GetPedestal();
265 nphot[i] += nratio*evtpix->GetNumPhotons();
266 perr[i] += Sqr(nratio*evtpix->GetErrorPhot());
267 pedrms[i] += Sqr(nratio*nppix.GetPedestalRms());
268
269 num++;
270 }
271
272 //
273 // Now the mean is calculated and the values rescaled back to the pixel area
274 //
275 nphot[i] /= num*ratio;
276 ped[i] /= num*ratio;
277 perr[i] = TMath::Sqrt(perr[i]/num)*ratio;
278 pedrms[i] = TMath::Sqrt(pedrms[i]/num)*ratio;
279
280 }
281
282 //
283 // Now the new pixel values are calculated and can be replaced in
284 // the corresponding containers
285 //
286 for (UShort_t i=0; i<entries; i++)
287 {
288 //
289 // Do not use blind neighbors
290 //
291 if (!fPixels->IsBlind(i))
292 continue;
293
294 //
295 // It must exist, we have created it in the loop before.
296 //
297 fEvt->GetPixById(i)->Set(nphot[i], perr[i]);
298 (*fPed)[i].Set(ped[i], pedrms[i]);
299 }
300
301 //
302 // Delete the intermediat arrays
303 //
304 delete nphot;
305 delete perr;
306 delete ped;
307 delete pedrms;
308}
309
310// --------------------------------------------------------------------------
311//
312// Removes all blind pixels from the analysis by setting their state
313// to unused.
314//
315void MBlindPixelCalc::Unmap() const
316{
317 const UShort_t entries = fEvt->GetNumPixels();
318
319 //
320 // remove the pixels in fPixelsIdx if they are set to be used,
321 // (set them to 'unused' state)
322 //
323 for (UShort_t i=0; i<entries; i++)
324 {
325 MCerPhotPix &pix = (*fEvt)[i];
326
327 if (fPixels->IsBlind(pix.GetPixId()))
328 pix.SetPixelUnused();
329 }
330}
331
332// --------------------------------------------------------------------------
333//
334// Treat the blind pixels
335//
336Int_t MBlindPixelCalc::Process()
337{
338 if (TESTBIT(fFlags, kUseInterpolation) && fGeomCam)
339 Interpolate();
340 else
341 Unmap();
342
343 return kTRUE;
344}
345
346// --------------------------------------------------------------------------
347//
348// - Check whether pixels to disable are available. If pixels are
349// given by the user nothing more is done.
350// - Otherwise try to determin the blind pixels from the starfield
351// given in MMcRunHeader.
352//
353Bool_t MBlindPixelCalc::ReInit(MParList *pList)
354{
355 if (TESTBIT(fFlags, kUseBlindPixels))
356 return kTRUE;
357
358 //
359 // If pixels are given by the user, we are already done
360 //
361 if (fPixelsIdx.GetSize() > 0)
362 return kTRUE;
363
364 //
365 // Delete the old array holding the blind pixels for the last file
366 //
367 fPixels->Clear();
368
369 if (!fGeomCam->InheritsFrom("MGeomCamMagic"))
370 {
371 *fLog << warn << "MBlindPixelCalc::ReInit: Warning - Starfield only implemented for Magic standard Camera... no action." << endl;
372 return kTRUE;
373 }
374
375 //
376 // Set as blind some particular pixels because of a particular
377 // Star Field of View.
378 //
379 MMcRunHeader *mcrun = (MMcRunHeader*)pList->FindObject("MMcRunHeader");
380 if (!mcrun)
381 {
382 *fLog << warn << "MBlindPixelCalc::ReInit: Warning - No run header available... no action." << endl;
383 return kTRUE;
384 }
385
386 Int_t rah, ram, ras;
387 Int_t ded, dem, des;
388 mcrun->GetStarFieldRa(&rah, &ram, &ras);
389 mcrun->GetStarFieldDec(&ded, &dem, &des);
390
391 if (rah!=5 || ram!=34 || ras!=32 || ded!=22 || dem!=0 || des!=55)
392 {
393 *fLog << warn << "Warning - Starfield unknown..." << endl;
394 return kTRUE;
395 }
396
397 //
398 // Case for Crab Nebula FOV
399 //
400 fPixels->SetPixelBlind(400);
401 fPixels->SetPixelBlind(401);
402 fPixels->SetPixelBlind(402);
403 fPixels->SetPixelBlind(437);
404 fPixels->SetPixelBlind(438);
405 fPixels->SetPixelBlind(439);
406
407 *fLog << inf;
408 *fLog << "FOV is centered at CRAB NEBULA: Setting 6 blind pixels" << endl;
409 *fLog << "to avoid bias values of analysis due to CRAB NEBULA:" << endl;
410 *fLog << " Pixels: 400, 401, 402, 437, 438, 439" << endl;
411
412 return kTRUE;
413}
414
415void MBlindPixelCalc::StreamPrimitive(ofstream &out) const
416{
417 out << " MBlindPixelCalc " << GetUniqueName();
418 if (fName!=gsDefName || fTitle!=gsDefTitle)
419 {
420 out << "(\"" << fName << "\"";
421 if (fTitle!=gsDefTitle)
422 out << ", \"" << fTitle << "\"";
423 out <<")";
424 }
425 out << ";" << endl;
426
427 if (TESTBIT(fFlags, kUseInterpolation))
428 out << " " << GetUniqueName() << ".SetUseInterpolation();" << endl;
429 if (TESTBIT(fFlags, kUseCentralPixel))
430 out << " " << GetUniqueName() << ".SetUseCentralPixel();" << endl;
431
432 if (fPixelsIdx.GetSize()==0)
433 return;
434
435 out << " {" << endl;
436 out << " TArrayS dummy;" << endl;
437 for (int i=0; i<fPixelsIdx.GetSize(); i++)
438 out << " dummy[" << i << "]=" << ((TArrayS)fPixelsIdx)[i] << ";" << endl;
439 out << " " << GetUniqueName() << ".SetPixels(dummy);" << endl;
440 out << " }" << endl;
441}
Note: See TracBrowser for help on using the repository browser.