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

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