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

Last change on this file since 4844 was 4819, checked in by tbretz, 20 years ago
*** empty log message ***
File size: 15.8 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(6), 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 if (TestFlags(~kCalcHillasSrc))
189 {
190 fCerPhotEvt = (MCerPhotEvt*)pList->FindObject(AddSerialNumber("MCerPhotEvt"));
191 if (!fCerPhotEvt)
192 {
193 *fLog << err << "MCerPhotEvt not found... aborting." << endl;
194 return kFALSE;
195 }
196 }
197
198 if (TestFlags(kCalcHillas|kCalcHillasExt|kCalcNewImagePar|kCalcConc))
199 {
200 fGeomCam = (MGeomCam*)pList->FindObject(AddSerialNumber("MGeomCam"));
201 if (!fGeomCam)
202 {
203 *fLog << err << "MGeomCam (Camera Geometry) missing in Parameter List... aborting." << endl;
204 return kFALSE;
205 }
206 }
207
208 // depend on whether MHillas is an in- or output container
209 if (TestFlag(kCalcHillas))
210 {
211 fHillas = (MHillas*)pList->FindCreateObj("MHillas", AddSerialNumber(fNameHillas));
212 if (!fHillas)
213 return kFALSE;
214 }
215
216 if (TestFlags(kCalcHillasExt|kCalcNewImagePar|kCalcConc|kCalcHillasSrc))
217 {
218 fHillas = (MHillas*)pList->FindObject(AddSerialNumber(fNameHillas), "MHillas");
219 if (!fHillas)
220 {
221 *fLog << err << fNameHillas << " [MHillas] not found... aborting." << endl;
222 return kFALSE;
223 }
224 }
225
226 // if enabled
227 if (TestFlag(kCalcHillasExt))
228 {
229 fHillasExt = (MHillasExt*)pList->FindCreateObj("MHillasExt", AddSerialNumber(fNameHillasExt));
230 if (!fHillasExt)
231 return kFALSE;
232 }
233
234 // if enabled
235 if (TestFlag(kCalcHillasSrc))
236 {
237 const MSrcPosCam *src = (MSrcPosCam*)pList->FindObject(AddSerialNumber(fNameSrcPosCam), "MSrcPosCam");
238 if (!src)
239 {
240 *fLog << warn << AddSerialNumber(fNameSrcPosCam) << " [MSrcPosCam] not found... creating default container." << endl;
241 src = (MSrcPosCam*)pList->FindCreateObj("MSrcPosCam", AddSerialNumber(fNameSrcPosCam));
242 }
243 if (!src)
244 return kFALSE;
245
246 fHillasSrc = (MHillasSrc*)pList->FindCreateObj("MHillasSrc", AddSerialNumber(fNameHillasSrc));
247 if (!fHillasSrc)
248 return kFALSE;
249
250 fHillasSrc->SetSrcPos(src);
251 }
252
253 // if enabled
254 if (TestFlag(kCalcNewImagePar))
255 {
256 fNewImgPar = (MNewImagePar*)pList->FindCreateObj("MNewImagePar", AddSerialNumber(fNameNewImagePar));
257 if (!fNewImgPar)
258 return kFALSE;
259 }
260
261 // if enabled
262 if (TestFlag(kCalcImagePar))
263 {
264 fImagePar = (MImagePar*)pList->FindCreateObj("MImagePar", AddSerialNumber(fNameImagePar));
265 if (!fImagePar)
266 return kFALSE;
267 }
268
269 // if enabled
270 if (TestFlag(kCalcConc))
271 {
272 fConc = (MConcentration*)pList->FindCreateObj("MConcentration", fNameConc);
273 if (!fConc)
274 return kFALSE;
275 }
276
277 fErrors.Reset();
278
279 Print();
280
281 return kTRUE;
282}
283
284// --------------------------------------------------------------------------
285//
286// If you want do complex descisions inside the calculations
287// we must move the calculation code inside this function
288//
289// If the calculation wasn't sucessfull skip this event
290//
291Int_t MHillasCalc::Process()
292{
293 if (TestFlag(kCalcHillas))
294 {
295 const Int_t rc = fHillas->Calc(*fGeomCam, *fCerPhotEvt, fIdxIsland);
296 if (rc<0 || rc>4)
297 {
298 *fLog << err << dbginf << "MHillas::Calc returned unknown error code!" << endl;
299 return kFALSE;
300 }
301 if (rc>0)
302 {
303 fErrors[rc]++;
304 return kCONTINUE;
305 }
306 }
307
308 if (TestFlag(kCalcHillasSrc))
309 {
310 const Int_t rc = fHillasSrc->Calc(*fHillas);
311 if (rc<0 || rc>2)
312 {
313 *fLog << err << dbginf << "MHillasSrc::Calc returned unknown error code!" << endl;
314 return kFALSE;
315 }
316 if (rc>0)
317 {
318 fErrors[rc+4]++;
319 return kCONTINUE;
320 }
321 }
322 fErrors[0]++;
323
324 if (TestFlag(kCalcHillasExt))
325 fHillasExt->Calc(*fGeomCam, *fCerPhotEvt, *fHillas, fIdxIsland);
326
327 if (TestFlag(kCalcImagePar))
328 fImagePar->Calc(*fCerPhotEvt);
329
330 if (TestFlag(kCalcNewImagePar))
331 fNewImgPar->Calc(*fGeomCam, *fCerPhotEvt, *fHillas, fIdxIsland);
332
333 if (TestFlag(kCalcConc))
334 fConc->Calc(*fGeomCam, *fCerPhotEvt, *fHillas);
335
336 return kTRUE;
337}
338
339// --------------------------------------------------------------------------
340//
341// This is used to print the output in the PostProcess. Later (having
342// millions of events) we may want to improve the output.
343//
344void MHillasCalc::PrintSkipped(int i, const char *str) const
345{
346 *fLog << " " << setw(7) << fErrors[i] << " (";
347 *fLog << setw(3) << (int)(100.*fErrors[i]/GetNumExecutions());
348 *fLog << "%) Evts skipped: " << str << endl;
349}
350
351// --------------------------------------------------------------------------
352//
353// Prints some statistics about the hillas calculation. The percentage
354// is calculated with respect to the number of executions of this task.
355//
356Int_t MHillasCalc::PostProcess()
357{
358 if (GetNumExecutions()==0)
359 return kTRUE;
360
361 if (!TestBit(kCalcHillas) && !TestBit(kCalcHillasSrc))
362 return kTRUE;
363
364 *fLog << inf << endl;
365 *fLog << GetDescriptor() << " execution statistics:" << endl;
366 *fLog << dec << setfill(' ');
367 if (TestBit(kCalcHillas))
368 {
369 PrintSkipped(1, "Less than 3 pixels (before cleaning)");
370 PrintSkipped(2, "Calculated Size == 0 (after cleaning)");
371 PrintSkipped(3, "Number of used pixels < 3");
372 PrintSkipped(4, "CorrXY==0");
373 }
374 if (TestBit(kCalcHillasSrc))
375 {
376 PrintSkipped(5, "Dist==0");
377 PrintSkipped(6, "Arg2==0");
378 }
379 *fLog << " " << (int)fErrors[0] << " (" << (int)(100.*fErrors[0]/GetNumExecutions()) << "%) Evts survived Hillas calculation!" << endl;
380 *fLog << endl;
381
382 return kTRUE;
383}
384
385// --------------------------------------------------------------------------
386//
387// Print 'Dataflow'
388//
389void MHillasCalc::Print(Option_t *o) const
390{
391 *fLog << inf << GetDescriptor() << " calculating:" << endl;
392 if (TestFlag(kCalcHillas))
393 *fLog << " - " << fNameHillas << " from MGeomCam, MCerPhotEvt and fIdxIsland=" << fIdxIsland << endl;
394 if (TestFlag(kCalcHillasSrc))
395 *fLog << " - " << fNameHillasSrc << " from " << fNameSrcPosCam << ", " << fNameHillas << " and fIdxIsland=" << fIdxIsland << endl;
396 if (TestFlag(kCalcHillasExt))
397 *fLog << " - " << fNameHillasExt << " from MGeomCam, MCerPhotEvt, " << fNameHillas << " and fIdxIsland=" << fIdxIsland << endl;
398 if (TestFlag(kCalcImagePar))
399 *fLog << " - " << fNameImagePar << " from MCerPhotEvt" << endl;
400 if (TestFlag(kCalcNewImagePar))
401 *fLog << " - " << fNameNewImagePar << " from MGeomCam, MCerPhotEvt, " << fNameHillas << " and fIdxIsland=" << fIdxIsland << endl;
402 if (TestFlag(kCalcConc))
403 *fLog << " - " << fNameConc << " from MGeomCam, MCerPhotEvt and " << fNameHillas << endl;
404}
405
406// --------------------------------------------------------------------------
407//
408// Implementation of SavePrimitive. Used to write the call to a constructor
409// to a macro. In the original root implementation it is used to write
410// gui elements to a macro-file.
411//
412void MHillasCalc::StreamPrimitive(ofstream &out) const
413{
414 out << " MHillasCalc " << GetUniqueName() << "(";
415 if (fName!=gsDefName || fTitle!=gsDefTitle)
416 {
417 out << ", \"" << fName << "\"";
418 if (fTitle!=gsDefTitle)
419 out << ", \"" << fTitle << "\"";
420 }
421 out << ");" << endl;
422
423 if (TestFlags(kCalcHillasExt|kCalcNewImagePar|kCalcConc|kCalcHillasSrc))
424 {
425 if (fNameHillas!=gsNameHillas)
426 out << " " << GetUniqueName() << ".SetNameHillas(\"" << fNameHillas << "\");" << endl;
427 }
428 if (TestFlag(kCalcHillasSrc) && fNameHillasSrc!=gsNameHillasSrc)
429 {
430 out << " " << GetUniqueName() << ".SetNameHillasSrc(\"" << fNameHillasSrc << "\");" << endl;
431 out << " " << GetUniqueName() << ".SetNameSrcPosCam(\"" << fNameSrcPosCam << "\");" << endl;
432 }
433 if (TestFlag(kCalcHillasExt) && fNameHillasExt!=gsNameHillasExt)
434 out << " " << GetUniqueName() << ".SetNameHillasExt(\"" << fNameHillasExt << "\");" << endl;
435 if (TestFlag(kCalcConc) && fNameConc!=gsNameConc)
436 out << " " << GetUniqueName() << ".SetNameConc(\"" << fNameConc << "\");" << endl;
437 if (TestFlag(kCalcImagePar) && fNameImagePar!=gsNameImagePar)
438 out << " " << GetUniqueName() << ".SetNameImagePar(\"" << fNameImagePar << "\");" << endl;
439 if (TestFlag(kCalcNewImagePar) && fNameNewImagePar!=gsNameNewImagePar)
440 out << " " << GetUniqueName() << ".SetNameNewImagePar(\"" << fNameNewImagePar << "\");" << endl;
441
442 if (!TestFlag(kCalcHillas))
443 out << " " << GetUniqueName() << ".Disable(MHilllasCalc::kCalcHillas);" << endl;
444 if (!TestFlag(kCalcHillasExt))
445 out << " " << GetUniqueName() << ".Disable(MHilllasCalc::kCalcHillasExt);" << endl;
446 if (!TestFlag(kCalcHillasSrc))
447 out << " " << GetUniqueName() << ".Disable(MHilllasCalc::kCalcHillasSrc);" << endl;
448 if (!TestFlag(kCalcNewImagePar))
449 out << " " << GetUniqueName() << ".Disable(MHilllasCalc::kCalcNewImagePar);" << endl;
450 if (!TestFlag(kCalcConc))
451 out << " " << GetUniqueName() << ".Disable(MHilllasCalc::kCalcConc);" << endl;
452 if (!TestFlag(kCalcImagePar))
453 out << " " << GetUniqueName() << ".Disable(MHilllasCalc::kCalcImagePar);" << endl;
454
455 if (fIdxIsland>=0)
456 out << " " << GetUniqueName() << ".SetNumIsland(" << fIdxIsland << ");" << endl;
457}
458
Note: See TracBrowser for help on using the repository browser.