source: trunk/Mars/mimage/MHillasCalc.cc

Last change on this file was 18751, checked in by tbretz, 8 years ago
Implemented the possibility to correct for the reflector's abberation - here the outward shift of the CoG of the reflected light of a point source.
File size: 17.0 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, 12/2000 <mailto:tbretz@astro.uni-wuerzburg.de>
19! Author(s): Harald Kornmayer, 1/2001
20!
21! Copyright: MAGIC Software Development, 2000-2003
22!
23!
24\* ======================================================================== */
25
26/////////////////////////////////////////////////////////////////////////////
27//
28// MHillasCalc
29// ===========
30//
31// This is a task to calculate the Hillas parameters from each event
32//
33//
34// Flags
35// --------
36//
37// By default all flags are set:
38//
39// To switch of the calculation you may use:
40// - Disable(MHillasCalc::kCalcHillas)
41// - Disable(MHillasCalc::kCalcHillasExt)
42// - Disable(MHillasCalc::kCalcNewImagePar)
43// - Disable(MHillasCalc::kCalcImagePar)
44// - Disable(MHillasCalc::kCalcImagePar2)
45// - Disable(MHillasCalc::kCalcSrcPosCam)
46// - Disable(MHillasCalc::kCalcConc)
47//
48// If the calculation of MHillas is switched off a container MHillas
49// in the parameter list is nevertheless necessary for the calculation
50// of some other containers, see below.
51//
52// If kCalcHillasSrc is set and no corresponding MSrcPosCam is found
53// in the parameter list an empty container (X=0, Y=0) is created.
54//
55//
56// Container names
57// -----------------
58//
59// The names of the containers to be used can be set with:
60// - SetNameHillas("NewName")
61// - SetNameHillasExt("NewName")
62// - SetNameNewImgPar("NewName")
63// - SetNameImagePar("NewName")
64// - SetNameImagePar2("NewName")
65// - SetNameSrcPosCam("NewName")
66// - SetNameConc("NewName")
67// - SetNameHillasSrc("NewName")
68//
69//
70// Islands
71// ---------
72//
73// You can change the islands for which the caluclations are done by:
74// - SetNumIsland()
75// The default is to use all used pixels (-1)
76//
77// fIdxIslands effects the calculations:
78// - kCalcHillas
79// - kCalcHillasExt
80// - kCalcNewImgPar
81// - kCalcNewImgPar2
82//
83//
84// Example
85// ---------
86//
87// MHillasCalc calc0; // calculate all image parameters except source dep.
88// MHillasCalc calc1; // calculate source dependant image parameters for 'Source'
89// MHillasCalc calc2; // calculate source dependant image parameters for 'AntiSource'
90// MHillasCalc calc3; // calculate hillas parameters only for biggest island
91// MHillasCalc calc4; // calculate hillas parameter for 2nd biggest island
92// // setup names of input-/output-containers
93// calc1.SetNameSrcPosCam("Source");
94// calc2.SetNameSrcPosCam("AntiSource");
95// calc1.SetNameHillasSrc("MHillasSource");
96// calc2.SetNameHillasSrc("MHillasAntiSource");
97// calc3.SetNameHillas("MHillas0");
98// calc4.SetNameHillas("MHillas1");
99// // setup calculations to be done
100// calc0.Disable(MHillasCalc::kCalcHillasSrc);
101// calc1.SetFlags(MHillasCalc::kCalcHillasSrc);
102// calc2.SetFlags(MHillasCalc::kCalcHillasSrc);
103// calc3.SetFlags(MHillasCalc::kCalcHillas);
104// calc4.SetFlags(MHillasCalc::kCalcHillas);
105// // choode index of island
106// calc3.SetNumIsland(0);
107// calc4.SetNumIsland(1);
108//
109// // setup tasklist
110// MTaskList list;
111// list.Add(&calc0);
112// list.Add(&calc1);
113// list.Add(&calc2);
114// list.Add(&calc3);
115// list.Add(&calc4);
116//
117//
118// Input/Output Containers
119// -------------------------
120//
121// 1) MGeomCam 5) MHillas 8) MImagePar
122// 2) MSignalCam 6) MHillasSrc 9) MNewImagePar
123// 3) MSrcPosCam 7) MHillasExt 10) MNewImagePar2
124// 4) fIdxIslands 11) MConcentration
125//
126// Flag | Input Container | Output
127// -----------------+-----------------+--------
128// kCalcHillas | 1 2 4 | 5
129// kCalcHillasSrc | 3 4 5 | 6
130// kCalcHillasExt | 1 2 4 5 | 7
131// kCalcImagePar | 2 | 8
132// kCalcNewImgPar | 1 2 4 5 | 9
133// kCalcNewImgPar2 | 1 2 4 | 10
134// kCalcConc | 1 2 5 | 11
135// -----------------+-----------------+--------
136//
137/////////////////////////////////////////////////////////////////////////////
138#include "MHillasCalc.h"
139
140#include <fstream> // StreamPrimitive
141
142#include "MParList.h"
143
144#include "MSignalCam.h"
145
146#include "MHillas.h"
147#include "MHillasExt.h"
148#include "MHillasSrc.h"
149#include "MImagePar.h"
150#include "MNewImagePar.h"
151#include "MNewImagePar2.h"
152#include "MConcentration.h"
153
154#include "MLog.h"
155#include "MLogManip.h"
156
157ClassImp(MHillasCalc);
158
159using namespace std;
160
161const TString MHillasCalc::gsDefName = "MHillasCalc";
162const TString MHillasCalc::gsDefTitle = "Calculate Hillas and other image parameters";
163
164const TString MHillasCalc::gsNameHillas = "MHillas"; // default name of the 'MHillas' container
165const TString MHillasCalc::gsNameHillasExt = "MHillasExt"; // default name of the 'MHillasExt' container
166const TString MHillasCalc::gsNameNewImagePar = "MNewImagePar"; // default name of the 'MNewImagePar' container
167const TString MHillasCalc::gsNameNewImagePar2= "MNewImagePar2"; // default name of the 'MNewImagePar2' container
168const TString MHillasCalc::gsNameConc = "MConcentration"; // default name of the 'MConcentration' container
169const TString MHillasCalc::gsNameImagePar = "MImagePar"; // default name of the 'MImagePar' container
170const TString MHillasCalc::gsNameHillasSrc = "MHillasSrc"; // default name of the 'MHillasSrc' container
171const TString MHillasCalc::gsNameSrcPosCam = "MSrcPosCam"; // default name of the 'MSrcPosCam' container
172
173// --------------------------------------------------------------------------
174//
175// Default constructor.
176//
177MHillasCalc::MHillasCalc(const char *name, const char *title)
178 : fNameHillas(gsNameHillas), fNameHillasExt(gsNameHillasExt),
179 fNameHillasSrc(gsNameHillasSrc), fNameSrcPosCam(gsNameSrcPosCam),
180 fNameConc(gsNameConc), fNameImagePar(gsNameImagePar),
181 fNameNewImagePar(gsNameNewImagePar),
182 fNameNewImagePar2(gsNameNewImagePar2),
183 fErrors(7), fFlags(0xff), fIdxIsland(-1), fAbberation(1)
184{
185 fName = name ? name : gsDefName.Data();
186 fTitle = title ? title : gsDefTitle.Data();
187}
188
189// --------------------------------------------------------------------------
190//
191// Check for in-/output containers, see class description
192//
193Int_t MHillasCalc::PreProcess(MParList *pList)
194{
195
196 if (TestFlags(~kCalcHillasSrc))
197 {
198 fCerPhotEvt = (MSignalCam*)pList->FindObject(AddSerialNumber("MSignalCam"));
199 if (!fCerPhotEvt)
200 {
201 *fLog << err << "MSignalCam not found... aborting." << endl;
202 return kFALSE;
203 }
204 }
205
206 if (TestFlags(kCalcHillas|kCalcHillasExt|kCalcNewImagePar|kCalcNewImagePar2|kCalcConc))
207 {
208 fGeomCam = (MGeomCam*)pList->FindObject(AddSerialNumber("MGeomCam"));
209 if (!fGeomCam)
210 {
211 *fLog << err << "MGeomCam (Camera Geometry) missing in Parameter List... aborting." << endl;
212 return kFALSE;
213 }
214 }
215
216 // depend on whether MHillas is an in- or output container
217 if (TestFlag(kCalcHillas))
218 {
219 fHillas = (MHillas*)pList->FindCreateObj("MHillas", AddSerialNumber(fNameHillas));
220 if (!fHillas)
221 return kFALSE;
222 }
223
224 if (TestFlags(kCalcHillasExt|kCalcNewImagePar|kCalcConc|kCalcHillasSrc))
225 {
226 fHillas = (MHillas*)pList->FindObject(AddSerialNumber(fNameHillas), "MHillas");
227 if (!fHillas)
228 {
229 *fLog << err << fNameHillas << " [MHillas] not found... aborting." << endl;
230 return kFALSE;
231 }
232 }
233
234 // if enabled
235 if (TestFlag(kCalcHillasExt))
236 {
237 fHillasExt = (MHillasExt*)pList->FindCreateObj("MHillasExt", AddSerialNumber(fNameHillasExt));
238 if (!fHillasExt)
239 return kFALSE;
240 }
241
242 // if enabled
243 if (TestFlag(kCalcHillasSrc))
244 {
245 const MSrcPosCam *src = (MSrcPosCam*)pList->FindObject(AddSerialNumber(fNameSrcPosCam), "MSrcPosCam");
246 if (!src)
247 {
248 *fLog << warn << AddSerialNumber(fNameSrcPosCam) << " [MSrcPosCam] not found... creating default container." << endl;
249 src = (MSrcPosCam*)pList->FindCreateObj("MSrcPosCam", AddSerialNumber(fNameSrcPosCam));
250 }
251 if (!src)
252 return kFALSE;
253
254 fHillasSrc = (MHillasSrc*)pList->FindCreateObj("MHillasSrc", AddSerialNumber(fNameHillasSrc));
255 if (!fHillasSrc)
256 return kFALSE;
257
258 fHillasSrc->SetSrcPos(src);
259 }
260
261 // if enabled
262 if (TestFlag(kCalcNewImagePar))
263 {
264 fNewImgPar = (MNewImagePar*)pList->FindCreateObj("MNewImagePar", AddSerialNumber(fNameNewImagePar));
265 if (!fNewImgPar)
266 return kFALSE;
267 }
268
269 // if enabled
270 if (TestFlag(kCalcNewImagePar2))
271 {
272 fNewImgPar2 = (MNewImagePar2*)pList->FindCreateObj("MNewImagePar2", AddSerialNumber(fNameNewImagePar2));
273 if (!fNewImgPar2)
274 return kFALSE;
275 }
276
277 // if enabled
278 if (TestFlag(kCalcImagePar))
279 {
280 fImagePar = (MImagePar*)pList->FindCreateObj("MImagePar", AddSerialNumber(fNameImagePar));
281 if (!fImagePar)
282 return kFALSE;
283 }
284
285 // if enabled
286 if (TestFlag(kCalcConc))
287 {
288 fConc = (MConcentration*)pList->FindCreateObj("MConcentration", fNameConc);
289 if (!fConc)
290 return kFALSE;
291 }
292
293 fErrors.Reset();
294
295 Print();
296
297 return kTRUE;
298}
299
300// --------------------------------------------------------------------------
301//
302// If you want do complex descisions inside the calculations
303// we must move the calculation code inside this function
304//
305// If the calculation wasn't sucessfull skip this event
306//
307Int_t MHillasCalc::Process()
308{
309 if (TestFlag(kCalcHillas))
310 {
311 const Int_t rc = fHillas->Calc(*fGeomCam, *fCerPhotEvt, fIdxIsland);
312 if (rc<0 || rc>4)
313 {
314 *fLog << err << dbginf << "MHillas::Calc returned unknown error code!" << endl;
315 return kFALSE;
316 }
317 if (rc>0)
318 {
319 fErrors[rc]++;
320 return kCONTINUE;
321 }
322 }
323
324 if (TestFlag(kCalcHillasSrc))
325 {
326 const Int_t rc = fHillasSrc->Calc(*fHillas, fAbberation);
327 if (rc<0 || rc>2)
328 {
329 *fLog << err << dbginf << "MHillasSrc::Calc returned unknown error code!" << endl;
330 return kFALSE;
331 }
332 if (rc>0)
333 {
334 fErrors[rc+4]++;
335 return kCONTINUE;
336 }
337 }
338 fErrors[0]++;
339
340 if (TestFlag(kCalcHillasExt))
341 fHillasExt->Calc(*fGeomCam, *fCerPhotEvt, *fHillas, fIdxIsland);
342
343 if (TestFlag(kCalcImagePar))
344 fImagePar->Calc(*fCerPhotEvt);
345
346 if (TestFlag(kCalcNewImagePar))
347 fNewImgPar->Calc(*fGeomCam, *fCerPhotEvt, *fHillas, fIdxIsland);
348
349 if (TestFlag(kCalcNewImagePar2))
350 fNewImgPar2->Calc(*fGeomCam, *fCerPhotEvt, fIdxIsland);
351
352 if (TestFlag(kCalcConc))
353 fConc->Calc(*fGeomCam, *fCerPhotEvt, *fHillas);
354
355 return kTRUE;
356}
357
358// --------------------------------------------------------------------------
359//
360// Prints some statistics about the hillas calculation. The percentage
361// is calculated with respect to the number of executions of this task.
362//
363Int_t MHillasCalc::PostProcess()
364{
365 if (GetNumExecutions()==0)
366 return kTRUE;
367
368 if (!TestFlag(kCalcHillas) && !TestFlag(kCalcHillasSrc))
369 return kTRUE;
370
371 *fLog << inf << endl;
372 *fLog << GetDescriptor() << " execution statistics:" << endl;
373 if (TestFlag(kCalcHillas))
374 {
375 PrintSkipped(fErrors[1], "0-Pixel Event (before cleaning, MC event?)");
376 PrintSkipped(fErrors[2], "Calculated Size == 0 (after cleaning)");
377 PrintSkipped(fErrors[3], "Number of used pixels < 3");
378 //PrintSkipped(fErrors[4], "CorrXY==0");
379 }
380 if (TestFlag(kCalcHillasSrc))
381 {
382 PrintSkipped(fErrors[5], "Dist==0");
383 PrintSkipped(fErrors[6], "Arg2==0");
384 }
385 *fLog << " " << (int)fErrors[0] << " (" << Form("%5.1f", 100.*fErrors[0]/GetNumExecutions()) << "%) Evts survived Hillas calculation!" << endl;
386 *fLog << endl;
387
388 return kTRUE;
389}
390
391// --------------------------------------------------------------------------
392//
393// Print 'Dataflow'
394//
395void MHillasCalc::Print(Option_t *o) const
396{
397 *fLog << inf << GetDescriptor() << " calculating:" << endl;
398 if (TestFlag(kCalcHillas))
399 *fLog << " - " << fNameHillas << " from MGeomCam, MSignalCam and fIdxIsland=" << fIdxIsland << endl;
400 if (TestFlag(kCalcHillasSrc))
401 *fLog << " - " << fNameHillasSrc << " from " << fNameSrcPosCam << ", " << fNameHillas << " and fIdxIsland=" << fIdxIsland << endl;
402 if (TestFlag(kCalcHillasExt))
403 *fLog << " - " << fNameHillasExt << " from MGeomCam, MSignalCam, " << fNameHillas << " and fIdxIsland=" << fIdxIsland << endl;
404 if (TestFlag(kCalcImagePar))
405 *fLog << " - " << fNameImagePar << " from MSignalCam" << endl;
406 if (TestFlag(kCalcNewImagePar))
407 *fLog << " - " << fNameNewImagePar << " from MGeomCam, MSignalCam, " << fNameHillas << " and fIdxIsland=" << fIdxIsland << endl;
408 if (TestFlag(kCalcNewImagePar2))
409 *fLog << " - " << fNameNewImagePar2 << " from MGeomCam, MSignalCam, " << fNameHillas << " and fIdxIsland=" << fIdxIsland << endl;
410 if (TestFlag(kCalcConc))
411 *fLog << " - " << fNameConc << " from MGeomCam, MSignalCam and " << fNameHillas << endl;
412}
413
414// --------------------------------------------------------------------------
415//
416// Implementation of SavePrimitive. Used to write the call to a constructor
417// to a macro. In the original root implementation it is used to write
418// gui elements to a macro-file.
419//
420void MHillasCalc::StreamPrimitive(ostream &out) const
421{
422 out << " MHillasCalc " << GetUniqueName() << "(";
423 if (fName!=gsDefName || fTitle!=gsDefTitle)
424 {
425 out << ", \"" << fName << "\"";
426 if (fTitle!=gsDefTitle)
427 out << ", \"" << fTitle << "\"";
428 }
429 out << ");" << endl;
430
431 if (TestFlags(kCalcHillasExt|kCalcNewImagePar|kCalcConc|kCalcHillasSrc))
432 {
433 if (fNameHillas!=gsNameHillas)
434 out << " " << GetUniqueName() << ".SetNameHillas(\"" << fNameHillas << "\");" << endl;
435 }
436 if (TestFlag(kCalcHillasSrc) && fNameHillasSrc!=gsNameHillasSrc)
437 {
438 out << " " << GetUniqueName() << ".SetNameHillasSrc(\"" << fNameHillasSrc << "\");" << endl;
439 out << " " << GetUniqueName() << ".SetNameSrcPosCam(\"" << fNameSrcPosCam << "\");" << endl;
440 }
441 if (TestFlag(kCalcHillasExt) && fNameHillasExt!=gsNameHillasExt)
442 out << " " << GetUniqueName() << ".SetNameHillasExt(\"" << fNameHillasExt << "\");" << endl;
443 if (TestFlag(kCalcConc) && fNameConc!=gsNameConc)
444 out << " " << GetUniqueName() << ".SetNameConc(\"" << fNameConc << "\");" << endl;
445 if (TestFlag(kCalcImagePar) && fNameImagePar!=gsNameImagePar)
446 out << " " << GetUniqueName() << ".SetNameImagePar(\"" << fNameImagePar << "\");" << endl;
447 if (TestFlag(kCalcNewImagePar) && fNameNewImagePar!=gsNameNewImagePar)
448 out << " " << GetUniqueName() << ".SetNameNewImagePar(\"" << fNameNewImagePar << "\");" << endl;
449 if (TestFlag(kCalcNewImagePar2) && fNameNewImagePar2!=gsNameNewImagePar2)
450 out << " " << GetUniqueName() << ".SetNameNewImagePar2(\"" << fNameNewImagePar2 << "\");" << endl;
451
452 if (!TestFlag(kCalcHillas))
453 out << " " << GetUniqueName() << ".Disable(MHilllasCalc::kCalcHillas);" << endl;
454 if (!TestFlag(kCalcHillasExt))
455 out << " " << GetUniqueName() << ".Disable(MHilllasCalc::kCalcHillasExt);" << endl;
456 if (!TestFlag(kCalcHillasSrc))
457 out << " " << GetUniqueName() << ".Disable(MHilllasCalc::kCalcHillasSrc);" << endl;
458 if (!TestFlag(kCalcNewImagePar))
459 out << " " << GetUniqueName() << ".Disable(MHilllasCalc::kCalcNewImagePar);" << endl;
460 if (!TestFlag(kCalcNewImagePar2))
461 out << " " << GetUniqueName() << ".Disable(MHilllasCalc::kCalcNewImagePar2);" << endl;
462 if (!TestFlag(kCalcConc))
463 out << " " << GetUniqueName() << ".Disable(MHilllasCalc::kCalcConc);" << endl;
464 if (!TestFlag(kCalcImagePar))
465 out << " " << GetUniqueName() << ".Disable(MHilllasCalc::kCalcImagePar);" << endl;
466
467 if (fIdxIsland>=0)
468 out << " " << GetUniqueName() << ".SetNumIsland(" << fIdxIsland << ");" << endl;
469}
470
471// --------------------------------------------------------------------------
472//
473// Read the setup from a TEnv, eg:
474// MHillasCalc.IdxIsland: -1
475//
476Int_t MHillasCalc::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
477{
478 Bool_t rc = kFALSE;
479 if (IsEnvDefined(env, prefix, "IdxIsland", print))
480 {
481 rc = kTRUE;
482 SetIdxIsland(GetEnvValue(env, prefix, "IdxIsland", fIdxIsland));
483 }
484 return rc;
485}
Note: See TracBrowser for help on using the repository browser.