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

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