source: trunk/MagicSoft/Mars/mimage/MHillasCalc.cc@ 6634

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