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

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