source: trunk/MagicSoft/Mars/mfilter/MFMagicCuts.cc@ 7189

Last change on this file since 7189 was 7173, checked in by tbretz, 19 years ago
*** empty log message ***
File size: 14.5 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, 03/2005 <mailto:tbretz@astro.uni-wuerzburg.de>
19!
20! Copyright: MAGIC Software Development, 2000-2005
21!
22!
23\* ======================================================================== */
24
25/////////////////////////////////////////////////////////////////////////////
26//
27// MFMagicCuts
28//
29/////////////////////////////////////////////////////////////////////////////
30#include "MFMagicCuts.h"
31
32#include <fstream>
33
34#include "MHMatrix.h"
35
36#include "MLog.h"
37#include "MLogManip.h"
38
39#include "MParList.h"
40
41#include "MGeomCam.h"
42#include "MHillasSrc.h"
43#include "MHillasExt.h"
44#include "MParameters.h"
45#include "MPointingPos.h"
46
47ClassImp(MFMagicCuts);
48
49using namespace std;
50
51// --------------------------------------------------------------------------
52//
53// constructor
54//
55MFMagicCuts::MFMagicCuts(const char *name, const char *title)
56 : fHil(0), fHilSrc(0), fHilAnti(0), fHilExt(0), fThetaSq(0),
57 fDisp(0), fPointing(0), fMatrix(0), fVariables(30), fThetaCut(kNone),
58 fHadronnessCut(kAll)
59{
60 fName = name ? name : "MFMagicCuts";
61 fTitle = title ? title : "Class to evaluate the MagicCuts";
62
63 fMatrix = NULL;
64
65 AddToBranchList("MHillas.fWidth");
66 AddToBranchList("MHillas.fLength");
67 AddToBranchList("MHillas.fSize");
68 AddToBranchList("MHillasSrc.fAlpha");
69 AddToBranchList("MHillasSrcAnti.fAlpha");
70 AddToBranchList("MPointingPos.fZd");
71 AddToBranchList("MHillasExt.fMaxDist");
72 AddToBranchList("MHillasExt.fM3Long");
73
74 fVariables[0] = 1.42547; // Xi
75 fVariables[1] = 0.233773; // Theta^2
76 fVariables[2] = 0.2539460; // Area[0]
77 fVariables[3] = 5.2149800; // Area[1]
78 fVariables[4] = 0.1139130; // Area[2]
79 fVariables[5] = -0.0889; // M3long
80 fVariables[6] = 0.18; // Xi(theta)
81}
82
83// --------------------------------------------------------------------------
84//
85// The theta cut value GetThetaSqCut() is propageted to the parameter list
86// as 'ThetaSquaredCut' as MParameterD so that it can be used somewhere else.
87//
88Int_t MFMagicCuts::PreProcess(MParList *pList)
89{
90 MGeomCam *cam = (MGeomCam*)pList->FindObject("MGeomCam");
91 if (!cam)
92 {
93 *fLog << err << "MGeomCam not found... aborting." << endl;
94 return kFALSE;
95 }
96
97 fMm2Deg = cam->GetConvMm2Deg();
98
99 fThetaSq = (MParameterD*)pList->FindCreateObj("MParameterD", "ThetaSquared");
100 if (!fThetaSq)
101 return kFALSE;
102 fDisp = (MParameterD*)pList->FindCreateObj("MParameterD", "Disp");
103 if (!fDisp)
104 return kFALSE;
105
106 // propagate Theta cut to the parameter list
107 // so that it can be used somewhere else
108 MParameterD *thetacut = (MParameterD*)pList->FindCreateObj("MParameterD", "ThetaSquaredCut");
109 if (!thetacut)
110 return kFALSE;
111 thetacut->SetVal(GetThetaSqCut());
112 //thetacut->SetReadyToSave();
113
114 MParameterD *m3lcut = (MParameterD*)pList->FindCreateObj("MParameterD", "M3lCut");
115 if (!m3lcut)
116 return kFALSE;
117 m3lcut->SetVal(fVariables[5]);
118
119 MParameterD *dispxi = (MParameterD*)pList->FindCreateObj("MParameterD", "DispXi");
120 if (!dispxi)
121 return kFALSE;
122 dispxi->SetVal(fVariables[0]);
123
124 dispxi = (MParameterD*)pList->FindCreateObj("MParameterD", "DispXiTheta");
125 if (!dispxi)
126 return kFALSE;
127 dispxi->SetVal(fVariables[6]);
128
129 Print();
130
131 if (fMatrix)
132 return kTRUE;
133
134 //-----------------------------------------------------------
135
136 fHil = (MHillas*)pList->FindObject("MHillas");
137 if (!fHil)
138 {
139 *fLog << err << "MHillas not found... aborting." << endl;
140 return kFALSE;
141 }
142
143 fHilExt = (MHillasExt*)pList->FindObject("MHillasExt");
144 if (!fHilExt)
145 {
146 *fLog << err << "MHillasExt not found... aborting." << endl;
147 return kFALSE;
148 }
149
150 fPointing = (MPointingPos*)pList->FindObject("MPointingPos");
151 if (!fPointing)
152 {
153 *fLog << err << "MPointingPos not found... aborting." << endl;
154 return kFALSE;
155 }
156
157 fHilSrc = (MHillasSrc*)pList->FindObject("MHillasSrc");
158 if (!fHilSrc)
159 {
160 *fLog << err << "MHillasSrc not found... aborting." << endl;
161 return kFALSE;
162 }
163
164 if (fThetaCut&kOff)
165 {
166 fHilAnti = (MHillasSrc*)pList->FindObject("MHillasSrcAnti", "MHillasSrc");
167 if (!fHilAnti)
168 {
169 *fLog << warn << "MHillasSrcAnti [MHillasSrc] not found... aborting." << endl;
170 return kFALSE;
171 }
172 }
173
174 return kTRUE;
175}
176
177// --------------------------------------------------------------------------
178//
179// Returns the mapped value from the Matrix
180//
181Double_t MFMagicCuts::GetVal(Int_t i) const
182{
183 return (*fMatrix)[fMap[i]];
184}
185
186// --------------------------------------------------------------------------
187//
188// You can use this function if you want to use a MHMatrix instead of the
189// given containers. This function adds all necessary columns to the
190// given matrix. Afterward you should fill the matrix with the corresponding
191// data (eg from a file by using MHMatrix::Fill). If you now loop
192// through the matrix (eg using MMatrixLoop) MEnergyEstParam::Process
193// will take the values from the matrix instead of the containers.
194//
195void MFMagicCuts::InitMapping(MHMatrix *mat)
196{
197 if (fMatrix)
198 return;
199
200 fMatrix = mat;
201
202 fMap[kESize] = fMatrix->AddColumn("log10(MHillas.fSize)");
203 fMap[kEArea] = fMatrix->AddColumn("MHillas.GetArea*MGeomCam.fConvMm2Deg*MGeomCam.fConvMm2Deg");
204 fMap[kEM3Long] = fMatrix->AddColumn("MHillasExt.fM3Long*sign(MHillasSrc.fCosDeltaAlpha)*MGeomCam.fConvMm2Deg");
205 fMap[kEWdivL] = fMatrix->AddColumn("MHillas.fWidth/MHillas.fLength");
206 fMap[kEZd] = fMatrix->AddColumn("MPointingPos.GetZdRad");
207
208 fMap[kEAlpha] = fMatrix->AddColumn("MHillasSrc.fAlpha");
209 fMap[kEDist] = fMatrix->AddColumn("MHillasSrc.fDist*MGeomCam.fConvMm2Deg");
210 //fMap[kDCA] = fMatrix->AddColumn("MHillasSrc.fDCA*MGeomCam.fConvMm2Deg");
211 if (fThetaCut&kOff)
212 {
213 fMap[kEAlphaAnti] = fMatrix->AddColumn("MHillasSrcAnti.fAlpha");
214 fMap[kEDistAnti] = fMatrix->AddColumn("MHillasSrcAnti.fDist*MGeomCam.fConvMm2Deg");
215 fMap[kEM3LongAnti] = fMatrix->AddColumn("MHillasExt.fM3Long*sign(MHillasSrcAnti.fCosDeltaAlpha)*MGeomCam.fConvMm2Deg");
216 //fMap[kDCAAnti] = fMatrix->AddColumn("MHillasSrcAnti.fDCA*MGeomCam.fConvMm2Deg");
217 }
218
219 //fMap[kLength] = fMatrix->AddColumn("MHillas.fLength*MGeomCam.fConvMm2Deg");
220 //fMap[kWidth] = fMatrix->AddColumn("MHillas.fWidth*MGeomCam.fConvMm2Deg");
221}
222
223// --------------------------------------------------------------------------
224//
225// Returns theta squared based on the formula:
226// p := c*(1-width/length)
227// return d*d + p*p - 2*d*p*cos(a*TMath::DegToRad());
228//
229// If par!=NULL p is stored in this MParameterD
230//
231Double_t MFMagicCuts::GetThetaSq(Double_t c, Double_t wl, Double_t d, Double_t a, MParameterD *par) const
232{
233 // wl = width/l [1]
234 // d = dist [deg]
235 // a = alpha [deg]
236 const Double_t p = c*(1-wl);
237 if (par)
238 par->SetVal(p);
239 return d*d + p*p - 2*d*p*cos(a*TMath::DegToRad());
240}
241
242// --------------------------------------------------------------------------
243//
244// Returns squared cut value in theta: fVariables[1]^2
245//
246Double_t MFMagicCuts::GetThetaSqCut() const
247{
248 return fVariables[1]*fVariables[1];
249}
250
251// ---------------------------------------------------------------------------
252//
253// Evaluate dynamical magic-cuts
254//
255Int_t MFMagicCuts::Process()
256{
257 // Get some variables
258 const Double_t wdivl = fMatrix ? GetVal(kEWdivL) : fHil->GetWidth()/fHil->GetLength();
259 const Double_t lgsize = fMatrix ? GetVal(kESize) : TMath::Log10(fHil->GetSize());
260 const Double_t zd = fMatrix ? GetVal(kEZd) : fPointing->GetZdRad();
261
262 // For simplicity
263 const Double_t *c = fVariables.GetArray();
264
265 // Value for Theta cut
266 const Double_t cut = GetThetaSqCut(); //c[1]*c[1];
267 const Double_t xi = (c[0]+c[8]*(lgsize-c[7])*(lgsize-c[7])) - c[6]*zd*zd;
268
269 // Default if we return before the end
270 fResult = kFALSE;
271
272 // ------------------- Most efficient cut -----------------------
273 // ---------------------- Theta cut ON --------------------------
274 const Double_t dist = fMatrix ? GetVal(kEDist) : fHilSrc->GetDist()*fMm2Deg;
275 const Double_t alpha = fMatrix ? GetVal(kEAlpha) : fHilSrc->GetAlpha();
276
277 const Double_t thetasq = GetThetaSq(xi, wdivl, dist, alpha, fDisp);
278
279 fThetaSq->SetVal(thetasq);
280
281 if (fThetaCut&kOn && thetasq >= cut)
282 return kTRUE;
283
284 // --------------------- Efficient cut -------------------------
285 // ---------------------- Area/M3l cut --------------------------
286 if (fHadronnessCut&kArea)
287 {
288 const Double_t area = fMatrix ? GetVal(kEArea) : fHil->GetArea()*fMm2Deg*fMm2Deg;
289 const Double_t A = lgsize>c[3] ? c[2] : c[2]*(1 - c[4]*(lgsize-c[3])*(lgsize-c[3]));
290 if (area>=A)
291 return kTRUE;
292 }
293 if (fHadronnessCut&kM3Long)
294 {
295 const Double_t m3long = fMatrix ? GetVal(kEM3Long) : fHilExt->GetM3Long()*TMath::Sign(fMm2Deg, fHilSrc->GetCosDeltaAlpha());
296 if (m3long<=c[5])
297 return kTRUE;
298 }
299
300 // ------------------- Least efficient cut -----------------------
301 // ---------------------- Theta cut OFF --------------------------
302 if (fThetaCut&kOff)
303 {
304 const Double_t distanti = fMatrix ? GetVal(kEDistAnti) : fHilAnti->GetDist()*fMm2Deg;
305 const Double_t alphaanti = fMatrix ? GetVal(kEAlphaAnti) : fHilAnti->GetAlpha();
306 const Double_t thetasqanti = GetThetaSq(xi, wdivl, distanti, alphaanti);
307 const Double_t m3lanti = fMatrix ? GetVal(kEM3LongAnti) : fHilExt->GetM3Long()*TMath::Sign(fMm2Deg, fHilAnti->GetCosDeltaAlpha());
308
309 if (thetasqanti<cut && (fHadronnessCut&kM3Long || m3lanti>c[5]))
310 return kTRUE;
311 }
312
313 fResult = kTRUE;
314
315 return kTRUE;
316}
317
318void MFMagicCuts::SetVariables(const TArrayD &arr)
319{
320 fVariables = arr;
321}
322
323TString MFMagicCuts::GetParam(Int_t i) const
324{
325 if (i>=fVariables.GetSize())
326 return "";
327
328 return Form("Param[%2d]=%+f", i, fVariables[i]);
329}
330
331void MFMagicCuts::Print(Option_t *o) const
332{
333 *fLog << all << GetDescriptor() << ":" << setiosflags(ios::left) << endl;
334
335 *fLog << "Using Theta cut: ";
336 switch (fThetaCut)
337 {
338 case kNone:
339 *fLog << "none" << endl;
340 break;
341 case kOn:
342 *fLog << "on" << endl;
343 break;
344 case kOff:
345 *fLog << "off" << endl;
346 break;
347 case kWobble:
348 *fLog << "on and off (wobble)" << endl;
349 break;
350 }
351 *fLog << "Using hadronness cut: ";
352 switch (fHadronnessCut)
353 {
354 case kNoCut:
355 *fLog << "none" << endl;
356 break;
357 case kArea:
358 *fLog << "area" << endl;
359 break;
360 case kM3Long:
361 *fLog << "m3long" << endl;
362 break;
363 case kAll:
364 *fLog << "all" << endl;
365 break;
366 }
367 *fLog << "Filter is" << (IsInverted()?"":" not") << " inverted." << endl;
368
369 const Int_t n = fVariables.GetSize();
370 for (int i=0; i<n; i+=3)
371 {
372 *fLog << setw(25) << GetParam(i);
373 *fLog << setw(25) << GetParam(i+1);
374 *fLog << setw(25) << GetParam(i+2);
375 *fLog << endl;
376 }
377 *fLog << setiosflags(ios::right);
378}
379
380Bool_t MFMagicCuts::CoefficentsRead(const char *fname)
381{
382 ifstream fin(fname);
383 if (!fin)
384 {
385 *fLog << err << "Cannot open file " << fname << ": ";
386 *fLog << strerror(errno) << endl;
387 return kFALSE;
388 }
389
390 for (int i=0; i<fVariables.GetSize(); i++)
391 {
392 fin >> fVariables[i];
393 if (!fin)
394 {
395 *fLog << err << "ERROR - Not enough coefficients in file " << fname << endl;
396 return kERROR;
397 }
398 }
399 return kTRUE;
400}
401
402Bool_t MFMagicCuts::CoefficentsWrite(const char *fname) const
403{
404 ofstream fout(fname);
405 if (!fout)
406 {
407 *fLog << err << "Cannot open file " << fname << ": ";
408 *fLog << strerror(errno) << endl;
409 return kFALSE;
410 }
411
412 for (int i=0; i<fVariables.GetSize(); i++)
413 fout << setw(11) << fVariables[i] << endl;
414
415 return kTRUE;
416}
417
418Int_t MFMagicCuts::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
419{
420 if (MFilter::ReadEnv(env, prefix, print)==kERROR)
421 return kERROR;
422
423
424 Bool_t rc = kFALSE;
425 if (IsEnvDefined(env, prefix, "ThetaCut", print))
426 {
427 TString str = GetEnvValue(env, prefix, "ThetaCut", "");
428 str.ToLower();
429 str=str.Strip(TString::kBoth);
430
431 if (str==(TString)"nocut" || str==(TString)"none")
432 fThetaCut = kNone;
433 if (str==(TString)"on")
434 fThetaCut = kOn;
435 if (str==(TString)"off")
436 fThetaCut = kOff;
437 if (str==(TString)"wobble")
438 fThetaCut = kWobble;
439 rc = kTRUE;
440 }
441 if (IsEnvDefined(env, prefix, "HadronnessCut", print))
442 {
443 TString str = GetEnvValue(env, prefix, "HadronnessCut", "");
444 str.ToLower();
445 str=str.Strip(TString::kBoth);
446
447 if (str==(TString)"nocut" || str==(TString)"none")
448 fHadronnessCut = kNoCut;
449 if (str==(TString)"area")
450 fHadronnessCut = kArea;
451 if (str==(TString)"m3long")
452 fHadronnessCut = kM3Long;
453 if (str==(TString)"all")
454 fHadronnessCut = kAll;
455
456 rc = kTRUE;
457 }
458
459 if (IsEnvDefined(env, prefix, "File", print))
460 {
461 const TString fname = GetEnvValue(env, prefix, "File", "");
462 if (!CoefficentsRead(fname))
463 return kERROR;
464
465 return kTRUE;
466 }
467
468 for (int i=0; i<fVariables.GetSize(); i++)
469 {
470 if (IsEnvDefined(env, prefix, Form("Param%d", i), print))
471 {
472 // It is important that the last argument is a floating point
473 fVariables[i] = GetEnvValue(env, prefix, Form("Param%d", i), 0.0);
474 rc = kTRUE;
475 }
476 }
477 return rc;
478 //return kTRUE; // means: can use default values
479 //return rc; // means: require something in resource file
480}
Note: See TracBrowser for help on using the repository browser.