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

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