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

Last change on this file since 4781 was 4711, checked in by tbretz, 20 years ago
*** empty log message ***
File size: 15.6 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(5), 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 memset(fErrors.GetArray(), 0, sizeof(Char_t)*fErrors.GetSize());
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 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 fErrors[rc]++;
302 if (rc>0)
303 return kCONTINUE;
304 }
305
306 if (TestFlag(kCalcHillasSrc))
307 {
308 if (!fHillasSrc->Calc(*fHillas))
309 {
310 fErrors[5]++;
311 return kCONTINUE;
312 }
313 }
314
315 if (TestFlag(kCalcHillasExt))
316 fHillasExt->Calc(*fGeomCam, *fCerPhotEvt, *fHillas, fIdxIsland);
317
318 if (TestFlag(kCalcImagePar))
319 fImagePar->Calc(*fCerPhotEvt);
320
321 if (TestFlag(kCalcNewImagePar))
322 fNewImgPar->Calc(*fGeomCam, *fCerPhotEvt, *fHillas, fIdxIsland);
323
324 if (TestFlag(kCalcConc))
325 fConc->Calc(*fGeomCam, *fCerPhotEvt, *fHillas);
326
327 return kTRUE;
328}
329
330// --------------------------------------------------------------------------
331//
332// This is used to print the output in the PostProcess. Later (having
333// millions of events) we may want to improve the output.
334//
335void MHillasCalc::PrintSkipped(int i, const char *str) const
336{
337 *fLog << " " << setw(7) << fErrors[i] << " (";
338 *fLog << setw(3) << (int)(100.*fErrors[i]/GetNumExecutions());
339 *fLog << "%) Evts skipped: " << str << endl;
340}
341
342// --------------------------------------------------------------------------
343//
344// Prints some statistics about the hillas calculation. The percentage
345// is calculated with respect to the number of executions of this task.
346//
347Int_t MHillasCalc::PostProcess()
348{
349 if (GetNumExecutions()==0)
350 return kTRUE;
351
352 if (!TestBit(kCalcHillas) && !TestBit(kCalcHillasSrc))
353 return kTRUE;
354
355 *fLog << inf << endl;
356 *fLog << GetDescriptor() << " execution statistics:" << endl;
357 *fLog << dec << setfill(' ');
358 if (TestBit(kCalcHillas))
359 {
360 PrintSkipped(1, "Less than 3 pixels (before cleaning)");
361 PrintSkipped(2, "Calculated Size == 0 (after cleaning)");
362 PrintSkipped(3, "Number of used pixels < 3");
363 PrintSkipped(4, "CorrXY==0");
364 }
365 if (TestBit(kCalcHillasSrc))
366 PrintSkipped(5, "Dist==0");
367 *fLog << " " << (int)fErrors[0] << " (" << (int)(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.